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INTRODUCTION 

This  program  builds  off  of  our  extensive  experience  in  using  electrical  properties  of  prostate  to  distinguish 
malignant  from  benign  tissues  [1-5]  and  specifically  stems  from  exciting  new  data  published  in  The  Prostate  [6] 
in  which  we  demonstrated  significant  electrical  property  differences  between  high-  and  low-grade  prostate 
cancer.  These  electrical  properties  are  influenced  by  a  tissue’s  intra-  and  extra-cellular  composition, 
morphology,  and  cellular  constituency,  and  we  have  hypothesize  that  it  is  possible  to  use  these  properties  to 
discriminate  between  normal,  low-grade,  and  high-grade  malignant  formations  in  a  clinical  setting.  While 
measuring  these  properties  by  direct  contact  with  the  tissues  is  possible  in  invasive  experiments,  it  is  desirable 
to  develop  methods  to  do  so  in  a  non-invasive  fashion.  To  date  our  group,  and  other  groups  around  the  world, 
have  investigated  Electrical  Impedance  Tomography  and  Microwave  Imaging,  two  techniques  which  are  limited 
in  resolution  by  the  underlying  physics.  The  primary  objective  of  this  current  program  is  to  develop  a  high- 
resolution  MR-based  approach  to  imaging  the  electrical  properties  of  prostate  with  the  intent  of  producing  a 
system  potentially  able  to  image  cancer  grade.  This  is  possible  by  leveraging  ultra-novel  developments  in  MRI, 
and  our  extensive  experience  in  developing  technologies  to  gauge  and  assess  the  utility  of  electrical  properties 
for  prostate  cancer  detection  [1-6]  and  assessment  and  in  developing  computer  algorithms  to  transform 
electromagnetic  data  into  electrical  property  images  of  the  prostate  [7-14].  Specifically,  we  are  attempting  to 
use  the  maps  of  MRI  RF  field  data  acquired  with  safe  and  fast  sequences  to  create  high-resolution  electrical 
property  images  of  the  prostate.  We  are  developing  this  novel  technology,  evaluating  it  in  an  ex  vivo  setting, 
and  finally  assessing  the  feasibility  of  employing  this  imaging  modality  in  a  routine  clinical  cohort  of  patients 
with  the  intent  of  having  a  significant  and  immediate  impact  on  clinical  practice.  By  developing  this  high- 
resolution  electrical  property  imaging  modality  we  expect  to  produce  highly  sensitive  and  specific  images  of 
cancer  grade  within  the  prostate  and  ultimately  better  guide  clinicians  in  distinguishing  aggressive  from 
indolent  disease. 

Much  of  the  second  year  of  this  program  has  focused  on  ex  vivo  prostate  imaging,  preliminary  ex  vivo  data 
analysis,  continued  developing  of  MR-EPT  conductivity  reconstructions,  and  preparing  for  in  vivo  data 
collections.  We  are  currently  in  a  No  Cost  Extension  year  in  which  we  will  focus  on  completing  ex  vivo  and  in 
vivo  data  collection  and  statistical  analysis  of  our  data.  In  addition  we  will  be  preparing  publications  and 
proposals  focused  more  heavily  on  clinical  data  acquisition  and  evaluation. 

BODY 

The  following  research  summary  is  presented  in  terms  of  the  approved  Statement  of  Work,  with  each  task 
being  discussed  separately.  When  appropriate,  detailed  discussion  is  referenced  to  manuscripts  published, 
submitted,  or  in  preparation  which  are  provided  in  the  Appendix.  Information  provided  in  our  previous  annual 
report  is  omitted  here  and  instead  we  note  the  tasks  that  are  completed  and  reference  our  2014  annual  report. 
Note  that  future  task  and  objectives  to  be  completed  are  marked  as  TBC. 

SPECIFIC  AIM  1:  TO  DEVELOP  MR-EPT  FOR  PROSTATE  IMAGING 
Major  Task  1 :  Develop  computation  toolbox  for  MR-EPT 

a)  Build  Matlab-based  toolbox  for  computing  electrical  property  images 

Completed  and  reported  on  in  our  2014  Annual  Report.  On-going  investigation  continue  in  order  to  further 
optimize  and  validate  the  approach  developed. 

We  have  fully  implemented  our  proposed  method  as  a  MATLAB  toolbox  for  MR-EPT  image  reconstruction 
as  described  in  our  2014  annual  report.  Over  the  past  year  we  have  conducted  additional  simulations  to 
explore  the  impact  measurement  noise  has  on  our  reconstructed  images.  An  example  analysis,  shown  in 
Figure  1,  demonstrates  that  the  inverse  approach  we  developed  performs  better  than  previously  developed 
direct  approaches.  This  analysis  is  further  described  in  Appendix  1. 
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Figure  1.  2D  cross-sections  of  three-dimensional  reconstructions  of  synthetic  phase  data  with  noise  levels  of  5%,  10%,  15%,  and  20%, 
associated  with  columns  1  to  4,  respectively.  Each  row  in  the  figure  represents  a  reconstruction  with  a  different  algorithm:  the  top  row 
shows  (for  different  levels  of  noise)  reconstructions  with  a  direct  approach  that  has  previously  been  developed.  The  reconstruction 
successfully  identifies  the  conductivity  distribution,  showing  a  left-to-right  transition.  The  second  and  third  row  of  the  figure  represent 
reconstructions  with  the  inverse  approach  we  have  developed  as  part  of  this  program,  using  Quadratic  Regularization  (QR)  and  Total 
Variation  (TV)  regularization,  respectively.  The  QR  algorithm  is  able  to  identify  the  left-to-right  change  in  conductivity  and  to  describe  it 
with  a  certain  degree  of  smoothness,  which  is  characteristic  of  this  type  of  regularization.  The  TV  algorithm  is  able  to  identify  and 
describe  with  pronounced  sharpness  the  left-to-right  change  in  conductivity.  In  this  specific  case  the  TV  algorithm  also  appears 
particularly  robust  to  noise.  The  inverse  based  algorithms  seem  to  fare  batter  in  the  presence  of  noise  compared  to  the  direct  algorithm, 
as  it  would  be  expected  from  the  fact  that  they  do  not  need  to  differentiate  input  phase  data  (noisy  data  ->  large  variation  in  derivatives). 
All  figure  are  represented  on  a  grayscale  ranging  from  0  to  2  S/m. 


b)  Build  Matlab-based  toolbox  for  specific  MR-based  field  of  views 

Completed  and  reported  on  in  our  2014  Annual  Report.  On-going  investigation/developments  continue  in 
order  to  further  optimize  and  validate  this  toolbox. 

We  have  developed  Matlab-based  functions  to  read  in  arbitrary  MRI  DICOM  and  .PAR/.REC  (Phillips 
format)  images  and  display  them  for  evaluation.  Additional  developments  have  included  the  ability  to  create 
multi-slice  images  of  the  MR  data  for  use  in  comparing  the  different  MR  variants  we  are  exploring.  This 
toolbox  was  briefly  described  in  our  2014  annual  report. 

Task  1  Milestones: 

1 .  Functional  toolbox  for  producing  MR-EPT  images  -  Completed 

Major  Task  2:  Optimize  MR-EPT  and  multi-parametric  MR  imaging  through  phantom  imaging 

a)  Optimize  MR-sequence  for  MR-EPT 

Completed  and  reported  on  in  our  2014  Annual  Report. 

b)  Perform  initial  tank-based  phantom  imaging  studies 

Completed  and  reported  on  in  our  2014  Annual  Report.  On-going  investigation  continue  in  order  to  further 
optimize  and  validate  the  approach  developed. 

A  number  of  additional  tank-based  phantom  studies  were  conducted  during  this  past  year  to  further 
demonstrate  that  our  MR-EPT  algorithms  are  able  to  accurately  estimate  the  internal  conductivity  of  a 
volume.  One  such  experiment  specifically  focused  on  validating  that  our  approach  is  able  to  accurately 
reconstruct  curved  surfaces  (some  of  our  previous  experiments  explored  linearly  varying  conductivity 
changes).  Figure  2  displays  a  curvilinear  gelatin  phantom  conducted  along  with  the  MR  magnitude  image 
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produced.  The  phase  data  recorded  from  this  experiment  was  reconstructed  using  both  QR  and  TV 
approaches  to  solve  the  inverse  problem  (Figure  3).  Both  approaches  accurately  reconstruct  the  curved 
surface  suggesting  that  these  techniques  extend  beyond  reconstructing  linearly  varying  conductivity 
distributions.  In  addition,  we  divide  the  phase  images  into  a  number  of  sub-domains  to  reduce  the 
computational  time  associated  with  image  reconstruction.  We  explored  this  and  the  time  associated  with 
using  different  numbers  of  sub-domains  (Figure  3  and  Table  1).  The  experimental  configuration  and 
analysis  are  further  described  in  Appendix  1. 


Phantom  Photograph 


MRI  Magnitude 


Spatial  Location  [Pixel  Number] 

Figure  2.  Curved  phantom  photographic  image  (left)  and  MRI  magnitude  image  (right).  This  phantom  was  built  with  a  similar  method 
used  to  create  our  previous  linear  phantoms,  by  creating  a  gelatin  slab  which  is  immersed  in  a  saline  solution.  The  gelatin  slab  has  a 
conductivity  of  of  1.8  S/m  and  Copper  sulfate  (CuS04)  was  added  to  it,  to  generate  contrast  in  the  MRI  magnitude  image.  The  saline 
solution  has  a  conductivity  of  4.1  S/m.  In  part  (A)  of  the  figure  is  indicated  a  white  edge  that  approximates  the  perimeter  of  a  region  of 
interest  which  was  selected  for  image  reconstruction.  Part  (B)  of  the  figure  shows  the  MRI  magnitude  corresponding  to  the  region  of 
interest  used  in  the  reconstructions  of  Figure  3. 


Spatial  Location 


Spatial  Location 


Spatial  Location 


Quadratic  Reconstruction  3x3 


Spatial  Location 


Figure  3.  2D  cross-sections  of  three-dimensional  reconstructions  of  the 
3  e  phantom  of  Figure  2.  This  figure  demonstrates  the  ability  of  the  inverse 
25 1  algorithms  to  reconstruct  curved  boundaries,  and  also  serves  to 
L  I  demonstrate  the  possibility  of  speeding  up  image  reconstruction  by 
5 1  splitting  the  image  domain  into  smaller  subdomains  to  be  reconstructed 
u  independently.  The  left  column  represents  reconstructions  with  the 
inverse  algorithm  with  Quadratic  Regularization  and  the  right  column 
reconstructions  using  Total  Variation  regularization.  Each  row  represents 
a  different  subdomain  splitting  for  image  reconstruction.  In  the  top  row 

MRI  slices  have  been  split  in  two  parts  vertically,  indicated  as  2  x  1 

3  e  splitting,  in  the  second  row  the  slices  have  been  split  in  two  parts 
2S§  vertically  and  horizontally,  indicated  as  2  x  2  splitting,  and  the  third  and 

L  I  fourth  represent  respectively  3x2  and  3x3  splittings.  Splitting  the 

5 1  original  image  domain  in  smaller  parts  allows  for  a  faster  reconstruction 
u  as  reported  in  Table  I.  Some  discontinuities  are  present  at  those  image 
locations  where  two  different  subdomains  meet.  This  results  from  the  fact 
0  that  different  subdomains  are  treated  as  separate  by  the  reconstruction 
and  no  continuity  is  enforced.  In  future  work  we  intend  to  introduce  some 
correlation  between  neighboring  subdomains,  through  regularization 
1 3  1  techniques,  which  should  reduce  or  eliminate  these  discontinuities.  In 
2.5 1  both  sets  of  reconstructed  images,  although  some  minor  discontinuities 
L  I  are  present  at  the  interface  between  subdomains  both  and  inverse  QR  and 
5 1  inverse  TV  algorithms  are  able  to  correctly  reproduce  the  curved 
u  interface.  All  figures  are  represented  on  a  grayscale  ranging  from  0  to  4 
“  S/m. 

Table  1.  Computational  Time  For  Different  Subdomain  Configurations 

I  s' 


Config. 

Time  QR  Recon  [s] 

Time  TV  Recon  [s] 

Jac.  Mem.  [MB] 

2x1 

50.9 

684.8 

162 

2x2 

8.8 

232.9 

39 

3x2 

3.5 

133.3 

16 

3x3 

1.5 

85.9 

7 

c)  Perform  anatomically  accurate  phantom  imaging  studies 
Completed  and  reported  on  in  our  2014  Annual  Report. 
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Task  2  Milestones: 

Note  that  while  the  below  tasks  are  complete,  we  expect  to  continue  conducting  phantom  experiments  over  the 
course  of  the  next  year  to  continue  to  improve  our  image  reconstruction  algorithms  and  better  understand  any 
image  artifacts  that  may  appear  in  our  clinical  data  acquisition. 

1 .  Validated  MR-EPT  algorithms  -  Completed 

2.  Fully  functional  protocol  for  obtaining  MR-based  images  in  a  single  serially  acquired  imaging  session  - 

Completed 

3.  1  peer-reviewed  publication  submitted  -  Completed,  see  Appendix 

Major  Task  3:  Submit  documents  for  IRB  and  MRMC  HRPO  approval 

a)  Draft  and  submit  IRB  protocol  revisions  and  new  protocol  submission 
Ex  vivo  Protocol:  Completed  during  last  annual  reporting  period. 

In  Vivo  Protocol:  We  submitted  a  protocol  modification  request  to  the  CCRC  in  July  2014,  requesting 
approval  to  conduct  the  in  vivo  portion  of  the  proposed  work  (to  take  place  starting  month  15  of  this 
program);  in  consultation  with  our  local  IRB,  we  felt  this  staged  approach  was  appropriate  for  this  protocol. 
This  was  approved  by  CCRC  in  August  2014.  We  submitted  a  protocol  modification  to  the  Committee  to 
Protect  Human  Subjects  (Dartmouth’s  IRB  of  record)  in  September  2014.  This  was  approved  in  October 
2014. 

b)  Draft  and  submit  documentation  for  MRMC  HRPO  approval 

Ex  vivo  Protocol:  Completed  during  last  annual  reporting  period 

In  Vivo  Protocol:  We  submitted  a  protocol  modification  request  to  MRMS  HRPO  in  October  of  2014.  This 
was  approved  in  November  of  2014. 

Task  3  Milestones: 

1.  Obtain  IRB  and  MRMC  HRPO  approval  for  ex  vivo  and  in  vivo  cohorts  -  ex  vivo  completed,  in  vivo 
completed 

SPECIFIC  AIM  2:  TO  EVALUATE  MR-EPT  IN  AN  EX  VIVO  COHORT  OF  PROSTATES 
Major  Task  1:  Optimize  ex  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  MR-EPT  and  multi-parametric  MR  sequences  of  ex  vivo  prostates 
This  task  has  been  complete  and  described  in  our  2014  Annual  Report. 

b)  Optimize  MR-EPT  sequences  and  algorithms  based  on  findings  in  this  initial  cohort 
This  task  has  been  complete  and  described  in  our  2014  Annual  Report. 

Task  1  Milestones: 

1.  Validation  that  our  MR-EPT  and  multi-parametric  MR  protocol  is  initially  optimized,  robust,  and 
repeatable  -  Completed 

Major  Task  2:  Evaluate  ex  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  multi-parametric  MR  sequences  of  ex  vivo  prostates 

Over  the  past  year  we  have  actively  recruited  patients  to  participate  in  our  ex  vivo  study.  We  have  been 
averaging  ~  2  cases  per  month.  To  date  we  have  imaged  30  ex  vivo  prostates  (20  more  ex  vivo  prostates 
will  be  imaged  over  the  course  of  the  next  year  to  complete  our  target  of  50  ex  vivo  prostates).  For  all 
cases  we  have  been  recording  B1  phase  and  magnitude  images  (for  MR-EPT  image  reconstruction),  T1 
spin  echo,  T1  turbo  spin  echo,  and  T2  turbo  spin  echo  images.  -  TBC 

b)  Perform  semi-quantitative  and  quantitative  analysis  of  ex  vivo  prostate  samples 

We  have  developed  Matlab  code  to  analyze  our  MR-EPT  conductivity  images,  MR  images  (other  variants), 
and  pathology  maps.  Specifically,  we  have  developed  techniques  to  extract  regions  of  interest  from  the 
MR-based  images  for  comparison  with  the  pathology  maps.  Figure  4  describes  the  approachdeveloped  to 
analyze  our  MR  images  in  comparison  to  the  pathology  maps  being  generated  to  identify  cancer 
distribution  within  the  prostate.  Mean  values  from  these  ROIs  will  be  used  as  our  metric  for  statistical 
analysis  of  the  data. 
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Figure  4.  Approach  developed  to  statistically  analyze  the  images  in  comparison  to  the  pathology  maps  generated  for  each  case.  An 
automated  segmentation  algorithm  was  developed  to  identify  the  prostate  boundary  in  the  Tlw  images  (regions  enclosed  in  blue  in  MR 
images).  The  ellipses  within  the  Path  Map  were  warped  using  a  Thin  Plate  Spline  (TPS)  algorithm  to  match  the  segmented  boundaries 
identified  in  the  Tlw  images.  The  TPS  algorithm  also  appropriately  deforms  the  regions  identified  as  cancer  within  the  Path  Maps  and 
overlays  these  islands  on  the  MR-based  images  (regions  enclosed  in  magenta  in  the  MR  images).  The  pixels  within  the  magenta  regions 
are  categorized  as  cancer  and  those  outside  the  magenta  regions,  but  within  the  blue  regions  are  categorized  as  benign  tissues.  The 
pixel  categorization  are  ultimately  used  to  statistically  compare  benign  and  cancer  regions. 


c)  Statistically  analyze  MR-based  images  and  pathological  metrics 

Data  generation  continues  to  be  in  process  (MR  images  &  Pathology  metrics).  Figure  5  shows  a  typical 
panel  of  images  we  create  for  each  prostate  imaged.  We  also  have  the  ability  to  produce  a  variety  of  MR- 
EPT  images  based  on  different  algorithms  we  developed  (and  described  in  previous  reports)  (see  Figure  6 
for  an  example).  We  have  fully  reconstructed  25  of  our  30  cases.  Upon  initial  qualitative  review  of  our 
images  (qualitative  evaluation  of  our  initial  19  cases)  we  observed  an  artifact  in  our  inverse  3d 
reconstruction  algorithm.  This  artifact  manifests  as  an  unexpected  conductivity  variation  in  the  axial 
direction  (see  for  example  Figure  7).  Based  on  our  initial  qualitative  assessment  we  have  decided  to 
include  in  our  initial  analysis  of  the  data  MR-EPT  images  reconstructed  using  the  Laplacian  algorithm 
because  this  variation  doesn’t  appear  to  be  as  strong  (Figure  8).  We  are  in  the  process  of  evaluating  why 
our  inverse-based  algorithm  is  generating  this  artifact,  but  believe  it  is  associated  with  the  lack  of  continuity 
maintained  within  our  sub-domain  approach  to  image  reconstruction.  An  interesting  observation  that  we 
observed  during  our  initial  qualitative  assessment  is  that  there  appears  to  be  a  different  conductivity  within 
the  peripheral  zone  of  the  prostate  (see  for  example  Figure  9).  This  peripheral  enhancement  was  clearly 
present  in  10  of  19  cases,  possibly  present  in  7  of  19  cases,  clearly  not  present  in  1  of  19  cases,  and  1 
case  was  not  included  due  to  image  artifacts.  Figures  1 0  and  1 1  show  a  single  slice  from  16  of  our  cases 
(note  that  the  Laplacian-based  reconstruction  algorithm  demonstrates  much  less  patient-to-patient 
variability). 
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Figure  5:  Typical  panel  of  data  images  recorded  and  to  be  used  for  statistical  analysis,  a  denotes  conductivity  images. 


Figure  6.  Typical  panel  of  conductivity  images  based  on  applying  different  MR-EPT  algorithms  we  have  developed.  This  particular  example  displays  a  single 
slice  from  a  patient.  Five  different  reconstruction  algorithms  (discussed  in  previous  reports)  were  used  to  produce  these  conductivity  images.  We  plan  to 
statistically  analyze  each  of  the  algorithms  when  evaluating  the  clinical  efficacy  of  MR-EPT. 
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Figure  7.  Panel  of  16  conductivity  images  acquired  for  16  of  our  ex  vivo  prostate  cases  (this  panel  shows  results  from  our  3D  inverse  approach  to  image 
reconstruction).  We  aim  to  continue  exploring  this  variability  in  magnitude  during  the  next  year. 
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Figure  8.  MR-EPT  conductivity  images  for  a  single  patient.  A)  pathology  map,  b)  estimated  conductivity  based  on  our  3d  inverse-based  reconstruction 
algorithm,  c)  estimated  conductivity  based  on  a  laplacian-based  reconstruction  algorithm.  Note  the  conductivity  variation  along  the  axial  direction  of  the 
prostate  (i.e.  there  appears  to  be  conductivity  variations,  from  high  to  low  conductivity  (light  to  dark),  in  each  of  the  panels)  for  the  inverse  reconstruction 
algorithm  based  images.  These  variations  are  not  present  in  the  Laplacian-based  solution. 


Example  Case  1  Example  Case  2 


Figure  9.  Example  cases  exhibiting  enhancement  in  the  prostate’s  peripheral  zone.  Red  arrow  points  to  the  feature  in  each  of  the  cases.  The  enhancement  is 
observed  along  the  axis  of  the  prostate  (i.e.  observed  in  multiple  slices). 


Figure  10.  Panel  of  16  conductivity  images  acquired  for  16  of  our  ex  vivo  prostate  cases  (this  panel  shows  results  from  our  3D  inverse  approach  to  image 
reconstru ction ) .  We  aim  to  explore  the  variability  in  magnitude  during  the  next  quarter  and  to  conduct  an  interim  analysis  compare  these  images  with  the 
pathology  maps  we  have  produced  for  each  case.  Not  the  significant  variation  in  the  mean  conductivity  values  in  each  of  these  cases.  Based  on  our  qualitative 
evaluation,  this  variation  is  likely  due  to  the  axial  artifact  described  in  Figure  3. 


Figure  11.  Panel  of  16  conductivity  images  acquired  for  16  of  our  ex  vivo  prostate  cases  (this  panel  shows  results  from  our  3D  Laplacian  approach  to  image 
reconstru  ction ) .  Note  that  the  patient  to  patient  variation  is  much  less  in  this  series  of  images  as  compared  to  Figure  5. 


Ijiterim  Statistical  Analysis 

We  have  performed  an  initial  interim  analysis  on  18  of  our  30  cases  specifically  exploring  any  differences 
between  benign  and  cancer  pixels  as  categorized  through  the  segmentation  and  TPS-base  deformation 
used  to  co-register  our  Pathology  Maps  with  our  MR  images.  We  have  explored  the  data  on  an  individual 
patient  basis  and  as  a  composite  group  using  both  T-tests  and  Kolmogorov-Smirnov  (KS)  test.  Significant 
differences  in  the  mean  pixel  values  between  benign  and  malignant  regions  were  identified  in  14  of  18,  13 
of  18,  13  of  18,  and  11  of  18  patients  when  the  magnitude,  phase,  conductivity  (Laplacian-based),  and 
conductivity  (inverse  approach)  images  were  analyzed,  respectively  (Figure  12).  In  the  majority  of  cases  (9 
of  13),  the  conductivity  of  cancer  was  found  to  be  significantly  less  than  that  of  benign  pixels  (when  the 
Laplacian-based  images  are  considered).  This  agrees  with  theory  and  previous  findings  at  lower 
frequencies  in  which  the  denser  cancer  tissues  impeded  current  flow  more  than  in  benign  tissues  (i.e.  have 
a  lower  conductivity).  When  data  from  all  images  were  combined,  statistically  significant  differences  were 
noted  for  all  images  types  (magnitude,  phase,  Laplacian-based  conductivity,  inverse-based  conductivity) 
(Figure  13).  Despite  the  significant  differences,  there  is  moderate  variability  in  the  data  presented.  In 
addition  to  continuing  to  conduct  this  analysis  as  more  data  is  accrued,  we  will  be  coupling  the  Gleason 
grade  to  each  of  the  patient  data  sets  in  order  to  explore  conductivity  variations  with  respect  to  tissue 
morphology  (i.e.  Gleason  Grade).  We  may  find  that  the  cases  in  which  conductivity  in  cancer  is  less  than 
that  in  benign  regions  is  correlated  to  the  Gleason  grade.  This  will  be  part  of  our  on-going  analysis. 
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Figure  12.  Comparison  of  mean  cancer  and  benign  pixels  for  each  patient.  Magnitude,  phase,  Laplacian-based  conductivity,  and  Inverse-based 
conductivity  images  are  evaluated.  T-tests  and  KS-tests  were  used  to  evaluate  if  the  means  between  the  two  tissues  types  were  significant.  Gray  shaded 
regions  denote  the  patients  in  which  the  differences  were  NOT  significant.  Significant  differences  were  found  for  all  remaining  unshaded  patients. 
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Figure  13.  Comparison  of  mean  cancer  and  benign  pixels  for  entire  cohort  of  patients  analyzed  (18  prostates).  Magnitude,  phase,  Laplacian-based 
conductivity,  and  Inverse-based  conductivity  images  are  evaluated.  T-tests  and  KS-tests  were  used  to  evaluate  if  the  means  between  the  two  tissues  types 
were  significant.  Significant  differences  were  found  for  all  MR-imaging  variants. 


Task  2  Milestones: 

1.  Assessment  of  the  clinical  potential  MR-EPT  combined  with  multi-parametric  MR  might  have  for 
prostate  imaging  -  in  process  (TBC) 

2.  1  peer-reviewed  publication  submitted  -  in  process  (TBC) 
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SPECIFIC  AIM  3:  TO  EVALUATE  MR-EPT IN  AN  IN  VIVO  COHORT  OF  PATIENTS 
Major  Task  1:  Evaluate  in  vivo  MR-EPT  and  multi-parametric  MR  imaging 

a)  Record  MR-EPT  and  multi-parametric  MR  sequences  of  in  vivo  and  ex  vivo  prostates 

Recruiting  patients  for  our  in  vivo  protocol  is  in  progress.  Over  the  last  year  we  have  imaged  4  volunteers 
without  prostate  cancer.  We  initially  recorded  our  MR-EPT  sequences  from  a  volunteer  without  cancer  to 
evaluate  the  images  prior  to  enrolling  prostate  cancer  patients.  In  our  initial  volunteer  images  we  observed 
significant  phase  noise  in  the  images.  We  originally  traced  the  noise  to  a  breach  in  the  shielded  room 
surrounding  the  MR  system  we  have  been  using  in  our  Advanced  Imaging  System.  The  breach  did  not 
adversely  effect  our  ex  vivo  study  since  the  imaging  volume  is  so  small.  When  a  person  lays  in  the  MR 
bore,  they  act  as  an  antenna  absorbing  the  noise  in  the  MR  room  resulting  in  image  artifacts  within  our 
phase  images.  A  new  breach  was  identified  during  a  subsequent  in  vivo  test,  conducted  in  February, 
causing  additional  image  artifacts.  This  breach  was  repaired  in  early  April.  We  were  able  to  record 
additional  data  from  a  human  volunteer  (without  prostate  cancer)  in  June  to  confirm  that  the  images  have 
sufficient  quality  for  MR-EPT  reconstruction  (Figure  14).  Note  that  the  phase  images  have  a  similar 
appearance  to  our  ex  vivo  images,  with  substantially  less  noise  than  we  observed  in  our  earlier  in  vivo 
images.  The  Laplacian-based  reconstructions  were  observed  to  be  more  smooth  than  the  inverse-based 
reconstructions.  We  will  continue  to  explore  the  low  spatial  frequency  noise  exhibited  in  the  inverse  images 
going  forward.  Despite  these  noisy  inverse-based  images,  we  believe  that  the  phase  images  are 
sufficiently  noise  free  for  us  to  initiate  enrolling  prostate  cancer  patients.  We  are  scheduled  to  conduct  our 
first  in  vivo  prostate  imaging  case  in  August.  We  plan  to  enroll  1  patient  in  August  and  then  2  patients  per 
month  from  September  until  March  in  order  to  meet  our  targeted  enrollment  of  15  patients. 
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Figure  14.  First  in  vivo  MR-EPT  images  of  human  prostate.  This  was  a  volunteer  without  cancer.  Laplacian  reconstruction  show  more  promise  than 
that  inverse  reconstruction  in  in  vivo  conditions.  The  noise  in  the  Inverse-based  approach  is  likely  due  to  the  sub-domain  used  to  reconstruct  the  images. 
The  MR  signals  are  lower  in  the  large  volumes  of  tissue  images  in  vivo  as  compared  to  the  ex  vivo  cases.  We  believe  that  the  Laplacian-based  approach 
may  be  sufficient  for  MR-EPT  imaging.  In  addition  to  collecting  data  from  prostate  cancer  patients,  we  will  continue  to  explore  this  imaging  modality  in 
benign  volunteers.  One  approach  we  are  going  to  explore  for  improving  signal  to  noise  ratio  is  to  put  a  water  bladder  around  the  waist  of  the  patient. 
The  water  will  provide  additional  signal  for  MR  phase  imaging  which  may  help  to  improve  overall  image  quality. 

b)  Perform  post-prostatectomy  pathological  assessment  of  extracted  prostates 
To  be  completed  during  the  next  of  this  program 

c)  Statistically  analyze  MR-based  images  and  pathological  metrics 
To  be  completed  during  the  next  year  of  this  program 

Task  1  Milestones: 

1.  Comparison  between  in  vivo  and  ex  vivo  MR-EPT  -  TBC 

2.  Preliminary  clinical  statistics  defining  utility  of  MR-EPT  combined  with  other  MR-imaging  variants  -  TBC 

3.  Initial  parameter  threshold  values  for  use  in  detecting  and  staging  prostate  cancer  -  TBC 

4.  1  peer-reviewed  publication  submitted  -  TBC 
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KEY  RESEARCH  ACCOMPLISHMENTS 

•  Assembled  a  database  of  MR  and  MR-EPT  images  from  30  ex  vivo  human  prostate  samples;  this  will 
ultimately  be  made  available  to  the  research  community  once  we  have  acquired  data  for  the  proposed 
50  ex  vivo  prostates. 

•  Developed  a  robust  approach  for  registering  pathological  maps  of  ex  vivo  prostate  section  to  MR  and 
MR-EPT  image  stacks 

•  Demonstrated  through  an  interim  statistical  analysis  of  MR-EPT  for  ex  vivo  prostate  imaging  that 
significant  differences  between  cancer  and  benign  regions  exist  within  the  conductivity  images 

•  Produced  the  first  ever  MR-EPT-based  conductivity  images  of  in  vivo  human  prostate 


REPORTABLE  OUTCOMES 
Manuscripts 

Borsic  A,  Perreard  I,  Mahara  A,  Halter  RJ,  “An  Inverse  Problems  Approach  to  MR-EPT  Image  Reconstruction,” 
IEEE  Transactions  on  Medical  Imaging,  accepted  with  minor  modifications  June  2015.  (Appendix  1  - 
accepted  manuscript) 

Perreard  I,  Borsic  A,  Mahara  A,  Halter  RJ,  “Towards  Magnetic  Resonance  -  Electrical  Properties  Tomography 
(MR-EPT)  for  Prostate  Imaging,”  Medical  Physics,  to  be  submitted  September  2015  (Appendix  2  -  initial 
outline/draft) 


CONCLUSION 

It  is  a  daily  challenge  for  clinicians  to  determine  whether  a  man  recently  diagnosed  with  prostate  cancer  has 
aggressive  disease  requiring  immediately  radical  therapy  or  indolent  disease  requiring  a  more  passive  watchful 
waiting  or  active  surveillance  approach.  This  program  is  focused  on  developing  Magnetic  Resonance  - 
Electrical  Property  Tomography  (MR-EPT)  specifically  for  prostate  imaging.  Over  the  past  year  we  have 
continued  to  optimize  our  MR-EPT  algorithms  for  estimating  the  prostate’s  electrical  conductivity  given 
magnetic  field  phase  and  magnitude  images  acquired  using  custom  MR  sequences.  We  have  taken  corrective 
action  to  ensure  that  our  MR  suite  is  sufficiently  shielded  to  reduce  phase  noise  in  our  images.  Additional 
simulations  and  phantom  experiments  have  been  conducted  to  explore  the  influence  of  noise  on  our  MR-EPT 
images  and  to  validate  that  we  can  image  objects  with  curvilinear  structures.  We  have  continued  to  enroll  men 
in  our  ex  vivo  study,  have  recorded  data  from  30  ex  vivo  prostates,  and  have  demonstrated  significant 
conductivity  difference  between  benign  and  cancerous  regions  based  on  our  MR-EPT  images.  Finally,  the  first 
ever  in  vivo  prostate  conductivity  images  were  generated  using  MR-EPT.  Over  the  course  of  the  next  year,  we 
will  be  focusing  primarily  on  clinical  data  acquisition  (both  ex  vivo  and  in  vivo  cohorts)  and  in  analyzing  images 
acquired  to  assess  the  potential  of  using  MR-EPT  (coupled  with  other  MR  imaging  variants)  to  distinguish 
between  aggressive  and  indolent  prostate  cancer. 
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An  Inverse  Problems  Approach  to  MR-EPT  Image 

Reconstruction 

A.  Borsic,  I.  Perreard,  A.  Mahara,  R.  J.  Halter 


Abstract — Magnetic  Resonance  -  Electrical  Properties  Tomog¬ 
raphy  (MR-EPT)  is  an  imaging  modality  that  maps  the  spatial 
distribution  of  the  electrical  conductivity  and  permittivity  within 
a  body  using  information  from  the  B\  RF  field  captured  in  a 
standard  clinical  MR  system.  The  presence  of  a  body  within  the 
scanner  alters  the  phase  and  amplitude  of  the  B\  field,  and  by 
mapping  these  alterations  it  is  possible  to  recover  the  electrical 
properties.  The  B\  field  is  time-harmonic,  and  its  interaction  with 
the  imaged  body  can  be  described  by  the  Helmholtz  equation. 
Approximations  to  this  equation  have  been  previously  used  to 
estimate  the  conductivity  and  permittivity  of  the  imaged  body  in 
terms  of  first  or  second  derivatives  of  B\  phase  and  B\  amplitude 
data,  respectively.  Using  these  same  approximations,  an  inverse 
approach  to  solving  the  MR-EPT  problem  is  presented  here  that 
leverages  a  forward  model  for  describing  the  magnitude  and 
phase  of  the  B\  field  within  the  imaging  domain,  and  a  fitting 
approach  for  estimating  the  electrical  properties  distribution.  The 
advantages  of  this  approach  are  that  1)  differentiation  of  the 
measured  data  is  not  required,  thus  reducing  noise  sensitivity, 
and  2)  different  regularization  schemes  can  be  adopted,  depend¬ 
ing  on  prior  knowledge  of  the  distribution  of  conductivity  or 
permittivity,  leading  to  improved  image  quality.  To  demonstrate 
the  developed  approach,  both  quadratic  and  Total  Variation  (TV) 
regularization  methods  were  implemented  and  evaluated  through 
numerical  simulation  and  experimentally  acquired  data.  The 
proposed  inverse  approach  to  MR-EPT  reconstruction  correctly 
identifies  contrasts  and  accurately  reconstructs  the  geometry  in 
both  simulations  and  experiments.  The  TV  regularized  scheme 
reconstructs  sharp  spatial  transitions,  which  are  difficult  to 
reconstruct  with  other,  more  traditional  approaches.  While  the 
computational  burden  of  this  method  is  higher  compared  to 
previous  approaches  when  properties  are  estimated  over  large 
imaging  domains,  a  parallel  computational  approach  is  possible 
by  subdividing  the  imaging  domain  into  several  subdomains.  The 
feasibility  of  this  approach  is  demonstrated  on  numerical  and 
experimental  data.  The  memory  requirements  and  computational 
burden  of  this  reconstruction  approach  are  related  to  the  size  of 
the  Jacobian  matrix,  which  is  proportional  to  the  sixth  power  of 
the  size  of  the  imaging  domain.  Sub-dividing  the  domain  into  a 
few  subdomains  leads  to  significant  gains  in  computational  speed. 

Keywords'.  Conductivity,  Permittivity,  Reconstruction,  Elec¬ 
trical  Properties  Tomography,  Magnetic  Resonance,  Inverse 
Problem,  Quadratic  Regularization,  Total  Variation  Regular¬ 
ization,  Primal  Dual  -  Interior  Point  Method 

I.  Introduction 

Magnetic  Resonance  -  Electrical  Properties  Tomography 
(MR-EPT)  is  an  imaging  modality  that  maps  the  spatial  dis¬ 
tribution  of  the  electrical  conductivity  and  permittivity  using 
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standard  clinical  MR  systems.  This  technique  exploits  B\  field 
perturbations  associated  with  the  objects  present  within  the 
bore  of  the  scanner.  These  field  perturbations  can  be  used  to 
estimate  the  spatial  distribution  of  electrical  properties  within 
the  objects  giving  rise  to  the  perturbations. 

Early  work  in  MR-EPT  was  published  in  1991  [1],  but  the 
technique  has  only  recently  seen  a  resurgence  of  exploration 
due  to  the  clinical  potential  that  the  electrical  properties 
may  offer  in  terms  of  tissue  contrast.  Several  groups  have 
focused  on  developing  novel  image  reconstruction  methods 
and  algorithms  for  this  purpose. 

Reconstruction  algorithms  can  be  classified  as  direct  meth¬ 
ods ,  from  which  measured  B i  information  is  directly  used  to 
map  the  electrical  properties  (EPs),  and  as  inverse  methods ,  in 
which  EPs  are  estimated  by  fitting  a  model  to  the  measured 
data.  The  initial  work  by  Haake  et.al.  [1]  is  based  on  using  the 
Helmholtz  equation  for  describing  the  time  harmonic  magnetic 
field  H  at  radio  frequencies  used  in  MRI  imaging: 


V2H  =  - 


^  x  (V  x  H) 
k 


(1) 


where  k  =  e  —  j(cr/uj),  and  uj  is  the  angular  frequency 
of  the  H  field,  cr  the  electrical  conductivity,  e  the  electri¬ 
cal  permittivity,  and  fi  the  magnetic  permeability.  Assuming 
that  a  and  e  are  piecewise  constant  or  slowly  varying  (i.e. 
V/c  ~  0),  the  second  term  on  the  right  hand  side  of  (1)  can 
be  neglected.  Considering  only  the  MRI-measurable  positive 
circularly  polarized  component  H+  of  the  RF  transmit  field, 
and  considering  k  as  isotropic,  one  obtains: 


-1  V2iT+ 

k~ 


(2) 


which  expresses  the  EPs  as  a  function  of  the  H+  field  and  its 
second  derivative.  One  method  for  measuring  H+  amplitude 
is  through  use  of  double-angle  mapping  techniques  [2],  H+ 
phase  is  assumed  to  be  half  of  the  spin-echo  phase,  an 
assumption  valid  when  one  transmit  and  one  receive  coil  are 
used,  and  the  sensitivity  patterns  of  the  two  coils  have  a  similar 
spatial  distribution  but  a  reverse  polarity  [3].  Typically,  in 
biological  tissues,  /i  is  considered  equal  to  the  permeability 
of  free  space  /io- 

Despite  the  simplifying  assumptions  (2),  this  approach  has 
successfully  been  used  to  fit  layered  models  [1]  and  to  produce 
electrical  properties  images  in  post-mortem  animals  [4] .  While 
this  approach  is  feasible,  the  2nd-order  differentiation  required 
for  computing  the  Laplacian  is  sensitive  to  noise,  and  therefore 
not  desirable.  In  2009,  Katscher  et.al.  [3]  used  the  Gauss 
theorem  to  propose  an  alternate  formulation  that  decreases  the 
differentiation  from  2nd  to  1st  order.  With  this  transformation, 
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the  conductivity  and  permittivity  are  proportional  to  surface 
integrals  of  the  gradient  of  H+  phase  and  amplitude,  over  the 
surface  of  an  arbitrary  block  of  pixels  where  the  electrical 
properties  are  assumed  to  be  approximately  constant.  This 
approach  has  been  successfully  demonstrated  in  in-vitro  and  in 
in-vivo  experiments  [3],  [5],  [6]  and  exhibits  decreased  noise 
sensitivity  compared  to  previously  considered  approaches. 

The  assumption  that  EPs  are  constant  or  slowly  varying  also 
does  not  broadly  apply  to  imaging  of  biological  tissues,  and 
the  errors  arising  from  violating  this  assumption  are  analyzed 
in  detail  by  Seo  et.al.  [7]. 

An  approach  based  on  measuring  two  sets  of  H+  data 
with  multi-channel  MRI  systems  was  proposed  by  Zhang  et.al. 
[8].  This  approach  reduces  artifacts  resulting  from  fast  spatial 
variations  of  EPs,  but  requires  measuring  components  of  H+ 
in  the  (x,y)  plane  of  the  scanner,  which  is  impractical  in 
clinical  applications,  where  only  the  2  component  along  the 
direction  of  the  static  magnetic  field  is  normally  available. 

Sodickson  et.al  [9]  derived  a  general  formulation  to  MR- 
EPT  that  does  not  make  assumptions  on  the  distribution  of 
EPs  and  that  is  suitable  for  multi-channel  systems.  In  this 
approach  conductivity,  permittivity,  phase,  phase  derivatives, 
magnetization,  and  magnetization  derivatives  are  assumed  to 
be  unknowns.  This  approach  has  produced  promising  numeri¬ 
cal  results,  however,  it  requires  computing  spatially-dependent 
second  derivatives  and  may  therefore  be  sensitive  to  noise 
in  experimental  applications.  An  approach  based  on  formu¬ 
lating  the  dependence  of  EPs  on  H+  via  the  convection- 
reaction  equation  has  been  developed  by  Hafalir  et.al  [10]. 
This  approach  does  not  make  any  restrictive  assumption  on 
the  distribution  of  EPs.  However,  this  method  also  requires 
first  and  second  derivatives  of  H+  be  computed  and  is 
therefore  potentially  sensitive  to  noise;  despite  this,  it  has 
been  demonstrated  numerically  and  experimentally  to  provide 
results  that  are  superior  to  standard  methods  based  on  (2). 

Inverse  approaches  to  MR-EPT  have  been  also  considered. 
An  algorithm  termed  CSI-EPT  (Contrast  Source  Inversion  - 
EPT)  has  been  proposed  by  Balidemaj  et.al  [11].  This  algo¬ 
rithm  does  not  make  any  assumption  regarding  the  distribution 
of  EPs  and  is  based  on  an  inverse  formulation,  where  the 
EPs  are  fit  to  the  data  using  a  contrast  source  approach.  The 
method  has  been  shown  in  numerical  simulations  to  produce 
accurate  and  detailed  reconstructions  of  a  pelvis  model.  Total 
Variation  [12]  regularization  has  also  been  adopted  in  order 
to  reduce  this  method’s  sensitivity  to  noise.  To  the  best  of 
our  knowledge  this  method  has  not  yet  been  demonstrated 
on  experimental  data.  A  model-based  based  approach  to 
conductivity  reconstruction  which  incorporates  regularization 
techniques  has  also  been  proposed  by  Ropella  et.al  [13].  This 
approach  is  based  on  inverting  the  Laplacian  in  (3)  in  the 
Fourier  domain,  and  has  demonstrated  better  image  quality 
compared  to  traditional  direct  approaches  in  phantoms  and  in 
in-vivo  data. 

In  this  manuscript  we  present  a  novel  MR-EPT  recon¬ 
struction  algorithm  resulting  from  further  development  and 
enhancement  of  initial  work  developed  by  the  authors  [14], 
[15].  The  algorithm  is  based  on  an  inverse  approach  applied 
to  (2).  Conductivity  and  permittivity  are  treated  separately  as 


in  [4],  and  taken  as  parameters  to  be  fitted.  A  forward  model 
is  developed  to  link  electrical  parameters  to  the  H+  data. 
A  Jacobian  matrix,  based  on  the  forward  model,  is  defined 
and  used  for  updating  the  model  parameters.  Regularization 
is  used  to  stabilize  the  inversion  for  parameter  estimation. 
This  inverse  approach  to  MR-EPT,  based  on  fitting  H+  data 
with  a  model,  has  the  general  advantage  compared  to  direct 
methods  based  on  (2),  and  with  respect  to  [9],  [10],  of 
not  requiring  differentiation  of  measured  H+  data,  and  is 
therefore  less  sensitive  to  noise.  The  approach  is  demonstrated 
to  successfully  reconstruct  noisy  numerical  and  experimental 
data.  Because  this  approach  casts  MR-EPT  reconstruction  in  a 
well  established  inverse  problem  framework,  two  well  known 
regularization  techniques,  Quadratic  Regularization  and  Total 
Variation  regularization  are  implemented  to  reduce  the  effects 
of  noise  and  to  stabilize  the  inversion.  Quadratic  Regular¬ 
ization  leads  to  smoother  reconstructed  images,  while  Total 
Variation  regularization  is  able  to  produce  sharper  images. 
This  method,  being  based  on  the  approximate  relationship 
(2),  which  might  result  in  artifacts  at  the  interface  of  highly 
enterogenous  boundaries  as  discussed  in  [7].  However,  this 
approach  is  common  in  MR-EPT,  and  has  been  demonstrated 
to  produce  meaningful  images  [3],  [5],  [6],  [16],  [17]. 

Algorithms  based  on  (2)  have  been  applied  to  breast  cancer 
detection,  showing  potential  of  MR-EPT  as  a  diagnostic  tool. 
In  this  context  prior  structural  information  available  from  T1 
and/or  T2  weighted  MRI  imaging  has  been  used  to  enhance  the 
reconstructed  EPs  by  weighting  differently  variations  along  the 
normal  and  tangential  direction  with  respect  to  the  expected 
prior  features  [16],  [17].  Our  approach  also  enables  incorporat¬ 
ing  prior  structural  information,  so  that  preferential  directions 
of  change  of  EP  can  be  embedded  into  the  regularization  (e.g. 
[18],  [19]).  In  general  we  believe  that  the  approach  developed 
here  has  advantages  over  algorithms  based  directly  on  (2)  since 
it  does  not  require  differentiation  of  H+ .  This  may  make  it 
more  suitable  for  clinical  applications  where  noise  is  present. 
It  is  difficult  to  offer  comparative  remarks  with  respect  to 
recent  inverse  algorithms  [9],  [10],  [11],  as  these  comparisons 
will  likely  require  application  of  the  different  algorithms  on 
particular  tests  cases. 

Reconstruction  of  EPs  with  the  inverse  approach  devel¬ 
oped  her  entails  considering  the  EPs  of  every  pixel  in  the 
image  as  an  unknown  to  be  estimated.  For  three-dimensional 
datasets  this  can  result  in  an  excessive  computational  burden. 
We  describe  methods  for  splitting  the  imaging  domain  that 
significantly  reduce  this  burden. 

In  Section  II  we  introduce  our  reconstruction  approach  and 
develop  a  forward  model.  In  Section  III  we  describe  an  inverse 
formulation  for  reconstructing  EPs,  and  in  Section  IV  we 
discuss  an  implementation  using  Quadratic  Regularization  and 
one  using  Total  Variation  regularization.  Section  V  reports 
numerical  experiments,  Section  VI  discusses  the  computa¬ 
tional  burden  of  the  methods  developed  and  offers  methods  for 
reducing  it,  Section  VII  reports  image  reconstruction  results 
from  physical  experiments.  Finally  concluding  remarks  are 
offered. 
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II.  Forward  Model 

The  MR-EPT  reconstruction  approach  developed  here  is 
based  on  an  inverse  formulation,  in  which  a  model  is  fit  to 
measured  data.  In  this  context,  the  relationship  linking  the  H+ 
field  data  to  the  electrical  conductivity,  a,  and  permittivity, 
e,  defines  the  forward  model.  These  relationships  can  be 
approximated  as  [4]: 


and 


<7  =  v24>(h+) 

Mo  u 

1  V2|i7+| 

\H+\ 


(3) 

(4) 


Equation  (3)  can  be  recognized  as  the  Poisson  equation: 


V20  =  fiQUJcr 


(5) 


where  for  simplicity  0  denotes  0(77 + ).  Equation  (4)  takes  the 
form  of  the  Helmholtz  equation: 

V2A  +  Mo<^2eA  =  0  (6) 

where  for  simplicity  A  denotes  |77+|. 

This  manuscript  focuses  on  reconstructing  conductivity 
from  (5).  Permittivity  can  be  reconstructed  from  (6)  with 
similar  methods,  which  are  not  developed  in  this  manuscript. 
Boundary  conditions  appropriate  to  the  physics  describing  the 
forward  model  need  to  be  specified  in  order  to  solve  (5).  This 
approach  ultimately  aims  to  match  a  measured  phase  0meas 
with  a  simulated  0  computed  from  (5);  Dirichlet  conditions 
defining  the  phase  on  the  domain  boundary  are  therefore  ideal 
for  this  case.  This  conditions  is  specified  as 


<t>{r)  =  (f>{r)meas  Vr  e  dQ,  (7) 


where  r  is  a  point  in  space  and  dQ  is  the  boundary  of  the 
imaging  domain.  0meas  is  perfectly  matched  at  the  boundary 
through  use  of  (7),  and  it  will  be  matched  point-by -point  inside 
the  domain  by  appropriate  adjustment  of  the  a  distribution. 
It  is  worth  noting  that  for  any  value  of  measured  data, 
4>meas>  the  condition  (7)  is  compatible  with  (5)  [20],  and 
therefore  a  solution  to  (5)  always  exists  and  is  unique  for 
those  boundary  conditions  [20].  Equations  (5)  and  (7)  define 
a  forward  model  for  MR-EPT  conductivity  reconstruction  via 
the  inverse  approach  developed  here. 


III.  Inverse  Formulation 

The  forward  model  linking  0  to  a  can  be  used  to  establish  an 
estimate  of  the  a  distribution  by  fitting  the  phase  0  predicted 
by  the  model  to  the  measured  phase  0meas.  As  an  example, 
this  fitting  can  be  optimized  in  the  least  squares  sense  as 

crrec  =  argmin||0((7)  -  </>meas  l|2  (8) 

where  arec  is  the  spatial  distribution  of  the  electrical  conduc¬ 
tivity  reconstructed  by  fitting  0(a)  to  0meas.  This  fitting  is 
accomplished  by  adjusting  a  to  minimize  the  discrepancy  be¬ 
tween  0(a)  and  0meas  in  the  L2-norm  sense;  the  dependence 
of  the  model  predicted  phase  0  on  a  is  explicitly  shown  in 
(8). 


To  account  for  noisy  phase  measurements,  a  regularization 
term  [21]  is  adopted  to  stabilize  the  inversion,  transforming 

(8)  into 

CTree  =  argmin  [\\<I>{(t)  -  <t>meas\\2  +  a^(cr)]  (9) 

where  a  is  a  scalar  Tikhonov  factor  controlling  the  amount  of 
regularization  and  4/ (a)  is  a  regularization  functional.  These 
functionals  are  often  quadratic  and  involve  first  or  second 
differential  operators  [21],  but  many  different  forms  beyond 
simple  differential  operators  have  been  proposed  in  literature. 
Two  different  functionals  often  used  in  inverse  problems, 
Quadratic  Regularization  and  Total  Variation  Regularization, 
are  described  and  implemented  here. 

IV.  Implementation 

In  order  to  implement  an  MR-EPT  reconstruction  based  on 

(9) ,  a  regularization  functional  needs  to  be  defined,  a  numeric 
method  for  computing  0(a)  as  defined  by  (5)  and  (7)  must  be 
implemented,  and  the  Jacobian  matrix  of  the  mapping  a  i-a 
0(a)  needs  to  be  computed. 

A.  Quadratic  Regularization 

In  a  first  implementation  example,  a  classical  quadratic 
functional  can  be  used  to  transform  (9)  into 

CTree  =  argmin  [\\<j)(<r)  -  <?WaS||2  +  all-MI2]  (10) 

where  a  is  discretized  on  the  same  pixel  grid  as  that  used 
by  the  MR  scanner  to  map  0meas.  This  converts  a  from  a 
continuous  function  into  a  finite  vector  of  discretized  values. 
The  matrix  L  is  a  regularization  matrix,  which  we  have  chosen 
to  be  the  Laplacian  of  the  conductivity  distribution,  a  relatively 
common  choice  [22],  [23]. 

By  applying  the  Newton-Raphson  method  to  (10),  an  update 
equation  for  the  conductivity  can  be  derived  as 


Sct  =  -  [JTJ  +  aLTL]  1 

[JT (cp(a)  -  (pmeas)  -  aLT L(a  -  a*)]  (11) 


where  J  is  the  Jacobian  matrix  of  the  mapping  a  i-a  0(a), 
and  a*  is  an  initial  conductivity  distribution.  A  uniform 
conductivity  can  be  used  as  an  initial  starting  distribution  a* 
and  updated  using  (11),  as  arec  =  a*  +  5cr.  A  single  update 
is  sufficient  in  this  case  since  the  forward  model  is  linear  in 
a.  As  a  result,  the  Newton-Raphson  method  finds  the  solution 
to  (10)  in  one  step.  In  order  to  apply  (11)  one  has  to  compute 
0(a)  and  the  Jacobian  matrix  J. 

The  forward  problem  a  ^a  0(a)  in  three  dimensions 
consists  of  solving 


<920  <920  <920 

d^  +  W+d^=^ 


(12) 


where  the  (x,y,x)  are  the  coordinates  in  the  axes  (x,y,z), 
which  are  defined  to  be  aligned  to  the  main  axes  of  the  image 
stack.  The  Partial  Differential  Equation  (12)  can  be  easily 
discretized  on  the  image  grid  using  Finite  Difference  schemes 
[24]  expressing  the  partial  derivatives  as 
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d24>  _  <t>{x  +  hx)  -  2 <j>{x)  +  <j>(x  -  hx) 

dx 2  h2 

d2</>  =  4>{y  +  hy)  -  2(/>(y)  +  0(j/  -  hy) 

dy2  h2 

d2(f)  4>(z  +  hz)  —  2<fi(z)  +  4>(z  —  hz) 

dz 2  h2z 


(13a) 

(13b) 

(13c) 


where  hx,hy,  hz  are  respectively  the  pixel  spacings  along 
the  x,y,z  axes  and  the  points  (x  —  hx),x,  (x  +  hx),  (y  — 
hy),y ,  (y  A  hy),  (z  —  hz),z ,  and  (z  +  zx)  translate  into  the 
indices  within  the  vector  of  discrete  phases  (j).  As  is  standard 
in  Finite  Differences,  the  partial  differential  equation  (12)  can 
be  transformed  into  a  linear  system 


A(j)  =  b 


(14) 


where  A  is  a  matrix  derived  from  applying  (13)  to  the  different 
pixels  in  the  image,  b  is  a  right  hand  side  (RHS)  derived  from 
the  evaluation  of  the  term  /ji ea  in  (12)  at  each  pixel  location 
in  the  image,  and  0  is  the  vector  of  computed  phase  values. 

Appropriate  boundary  conditions  need  to  be  applied  to  the 
linear  system  in  (14).  The  following  approach  is  used  here 

•  Boundary  Pixels:  Dirichlet  boundary  conditions  are  es¬ 
tablished  for  all  the  pixels  on  the  boundary  surface  of 
the  image  stack  by  placing  a  1  on  the  diagonal  element 

where  i  is  the  index  of  each  of  these  pixels,  and 
setting  b(i)  =  4>{i)meaa- 

•  Internal  Pixels:  for  all  the  pixels  internal  to  the  do¬ 
main  (defined  as  being  at  least  1  pixel  away  from  the 
boundary),  (13)  is  applied  and  the  RHS  is  computed  as 
b(i)  =  /jLja(i),  where  i  is  the  index  of  the  considered 
pixel. 

Solving  the  linear  system  (14)  results  in  the  vector  <j)  of 
computed  phase  values. 

A  similar  approach  is  used  for  building  the  matrix  L,  which 
is  a  Laplacian  operator.  The  purpose  of  L  is  to  limit  high 
frequency  spatial  variations  in  the  image,  which  typically  arise 
from  noise  in  the  data.  The  term  ||Lcr||2  takes  large  values  for 
large  spatial  variations  or  for  high  frequency  spatial  variations 
in  cr,  which  penalizes  their  presence  in  the  reconstructed  image 
(10).  Elements  of  L  corresponding  to  the  interior  domain  are 
computed  using  the  same  finite  differences  applied  to  solve 
the  forward  problem  (13).  Elements  of  L  corresponding  to 
the  domain  boundary  are  defined  using  mirroring  boundary 
conditions.  Specifically,  if  the  term  (x  +  hx)  in  (13a)  falls 
outside  the  domain,  (x  +  hx)  is  assumed  equal  to  (x  —  hx) 
and  the  second  derivative  for  that  image  location  is  expressed 
as: 

d2(j>  =  —2<(>(x)  +  2 4>{x  -  hx) 

dx 2  *  h%  '  (  ’ 

Similar  conditions  are  used  for  derivatives  in  y  and  z. 

The  Jacobian  matrix  in  (11)  can  be  computed  from  (14) 
using  the  following  matrix  identities 


d(A  16)  |  db 

da  =  ~do 


(16) 


where  £  is  the  derivative  of  b  with  respect  to  cr.  In  this  case, 


£(i)  =0  for  all  indices  i  corresponding  to  boundary  pixels  (for 
which  the  Dirichlet  condition  has  been  set)  and  £(z)  =  yuj  for 
all  the  indices  corresponding  to  interior  pixels  where  b(i)  is  set 
to  b(i)  =  fjLJcr(i).  Equations  (14)  and  (16)  enable  calculation 
of  0  and  J,  which  are  in  turn  used  in  (11)  to  compute  arec. 

As  shown  later,  the  above  procedure  results  in  successful 
reconstructions  on  synthetic  and  experimental  data.  The  use  of 
the  quadratic  regularization  functional  produces  conductivity 
profiles  that  are  relatively  smooth.  The  benefit  of  the  inverse 
formulation  developed  here  (9)  is  that  different  functional 
terms  can  be  used  for  regularization.  In  the  next  subsection, 
Total  Variation  is  introduced  as  an  alternative  regularization 
term. 

B.  Total  Variation  Regularization 

In  the  previous  sections,  a  framework  for  MR-EPT  image 
reconstruction  based  on  inverse  problems  was  developed. 
One  benefit  of  this  framework  is  that  different  regularization 
terms  can  be  chosen.  Regularization  functionals  affect  how 
the  reconstructed  image  is  smoothed  and  different  choices 
are  appropriate  for  different  situations.  Total  Variation  (TV) 
is  a  relatively  novel  form  of  regularization  that  results  in 
images  with  sharper  conductivity  transitions  as  compared  to 
Quadratic  Regularization  approaches  (e.g.  11).  Reconstructing 
sharp  image  transitions  can  be  challenging  in  MR-EPT.  In 
direct  approaches  derivatives  are  estimated  on  multiple  points 
(e.g.  5,7,  or  9)  to  reduce  noise  sensitivity,  but  this  results  in 
smoothing.  In  the  inverse  approach  described  above  the  regu¬ 
larization  smooths  the  reconstructed  image,  again  for  reducing 
sensitivity  to  noise.  The  use  of  TV  regularization  instead 
allows  reconstructing  fast  spatial  variations  more  accurately. 

The  TV-based  inverse  formulation  for  (9)  is  expressed  as 

arec  =  argmin  [||0(cr)  -  4>meas\\2  +  «  TV (a)]  (17) 

where  the  TV  functional  is  defined  as  TV  (a)  = 

|V(<r)|df2  and  Q  is  the  imaging  domain.  The  sharper  re¬ 
constructions  possible  with  TV  regularization  arise  because  the 
TV  functional  remains  finite  for  step  changes,  while  quadratic 
functionals  like  fQ  |V(cr)|2<iD,  or,  \\72(a)\2dQ  (common 
quadratic  regularization  functionals)  tend  towards  infinity.  As 
a  result  of  the  large  quadratic  functional  values  associated 
with  fast  spatial  changes  in  conductivity,  these  profiles  (i.e. 
step  changes  in  conductivity)  are  penalized,  rendering  the 
reconstructions  smoother.  More  detailed  discussion  of  the  TV 
function  properties  are  presented  in  [12]. 

While  the  use  of  TV  is  desirable  for  reconstructing  sharp 
variations,  the  image  reconstruction  expressed  by  (17)  is  a  non- 
differentiable  optimization  problem,  and  special  techniques 
need  to  be  employed  to  minimize  \\(j)(a)  —  (frmeas  II 2 TV  (a) 
when  acting  on  a.  In  this  case,  a  Primal  Dual  -  Interior 
Point  framework  developed  in  [25],  [26]  for  optimizing  (17) 
was  used.  Specifically,  the  algorithm  named  “PD-IPM  -  L2- 
L1  Norm”  reported  in  [25]  is  used.  The  MR-EPT  forward 
model  (j){a )  and  Jacobian  matrix  J  developed  in  section  IV- A 
are  input  into  the  optimization  algorithm,  resulting  in  arec  as 
defined  in  (17). 

In  Sections  V  and  VII,  numerical  and  experimental  data 
are  used  respectively  to  compare  these  two  regularization 
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approaches  and  demonstrate  how  the  smoothing  /  sharpening 
characteristics  of  the  reconstruction  can  be  tuned  by  the 
regularization  functional. 

V.  Numerical  Experiments 

A  number  of  numerical  experiments  were  conducted  to 
validate  this  inverse  formulation  for  MR-EPT  conductivity 
reconstruction.  The  simulated  volume  of  interest  consists  of 
a  cubic  block  of  20x20x20  millimeters  split  in  half  by 
passing  a  plane  through  the  center  of  the  cube.  One  half  of 
the  cube  is  set  to  have  a  conductivity  of  1  Sm-1  while  the 
other  half  is  set  to  2  Sm_1,  as  shown  in  Figure  1.  While 
these  simulated  conductivities  are  generally  larger  than  typical 
physiological  values,  they  serve  as  a  good  test  do  demonstrate 
that  the  algorithms  can  reconstruct  a  unit  step  change  from 
1  to  2  Sm-1  at  the  interface  between  the  two  halves  of  the 
numerical  phantom. 

The  simulated  cube  was  discretized  into  lxlxl  mm  volume 
elements,  and  Finite  Differences  were  used  to  compute  the 
MRI  phase  using  equations  (12)  to  (14),  and  assuming  a  3 
Tesla  static  magnet  strength.  A  Dirichlet  boundary  condition 
of  (j){r)  =  0  Vr  £  dCl  was  assumed  on  the  boundary.  As 
discussed  in  Section  II,  Dirichlet  boundary  conditions  can  be 
used  for  matching  a  measured  phase  at  the  domain  boundary 
while  computing  the  forward  solution.  In  the  absence  of 
boundary  data,  the  condition  </>(r)  =  0  Vr  G  dCl  is  a  practical 
condition  that  is  compatible  with  the  Poisson  equation  and 
results  in  a  unique  solution.  This  choice  does  not  alter  the 
reconstructed  conductivity,  which  depends  on  the  Laplacian  of 
the  phase  -  the  value  of  which  is  enforced  within  the  domain 
by  the  Poisson  equation  itself. 

Figure  2  (A)  shows  a  2D  cross  section  of  the  computed 
synthetic  phase;  Figure  2  (B)  shows  the  same  synthetic  phase 
with  20%  additive  Gaussian  noise,  as  discussed  later  and  as 
used  in  the  reconstructions.  The  curvature  (i.e.  Faplacian)  of 
the  phase  is  more  pronounced  in  the  left  part  of  the  domain 
due  to  the  higher  conductivity  (2Sm-1)  within  this  region. 

The  simulated  phase  data  was  input  into  three  different 
reconstruction  algorithms.  The  first  is  an  implementation  of 
the  direct  approach  developed  by  Katscher  et.  al.  [3].  This 
algorithm  is  based  on  numerically  solving  (2)  and  using  a 
volume  of  integration  to  improve  robustness  to  noise.  Specif¬ 
ically,  the  order  of  differentiation  is  decreased  from  2nd  to 
1st  order  using  the  Gauss  theorem  to  convert  the  volume 
integral  of  (2)  to  a  surface  integral  which  has  only  first 
derivatives  of  the  field  variables.  This  represents  an  approach 
that  directly  computes  the  output  conductivity  as  a  function  of 
the  input  phase  distribution  by  taking  1st  order  derivatives  and 
computing  their  integral  over  specific  integration  volumes.  In 
our  implementation  of  this  algorithm,  an  integration  volume 
of  5  x  5  x  5  pixels  was  defined  and  phase  derivatives  were 
estimated  using  Savitzky-Golay  [27]  filters  involving  7  points 
(three  points  per  side  of  the  considered  pixel).  Katscher  et. 
al.  uses  similar  integration  volumes,  and  5  to  9  points  for 
derivative  estimation.  We  found  that,  for  our  numerical  data, 
using  only  7  points  represents  a  good  compromise  between 
image  sharpness  and  noise  sensitivity,  while  larger  numbers  of 


points  lead  to  smoother  images  with  limited  ability  to  identify 
sharp  transitions.  The  second  and  third  algorithms  are  imple¬ 
mentations  of  the  inverse  formulation  approach  developed  in 
this  manuscript  with  Quadratic  Regularization  (10)  and  with 
TV  regularization  (17),  respectively. 

Reconstructions  for  the  three  different  algorithms  are  shown 
in  Figure  3  using  a  fixed  gray-scale  for  all  figures  spanning 
conductivity  values  from  0.5  to  2  Sm-1.  The  top  row  of  the 
figure  represents  reconstruction  with  our  own  implementation 
of  the  algorithm  proposed  in  [3],  the  second  row  represents 
reconstruction  with  the  inverse  algorithm  using  Quadratic 
Regularization,  and  the  last  row  reconstructions  with  the 
inverse  algorithm  using  TV  regularization.  Different  levels  of 
simulated  noise  were  added  to  the  phase  data  for  each  column 
of  the  reconstructions.  Specifically,  noise  levels  of  5%,  10%, 
15%,  and  20%  were  added  to  columns  1  to  4,  respectively. 
Noise  was  generated  by  extracting  values  from  a  Gaussian 
random  distribution.  These  values  were  scaled  to  generate  a 
particular  percent  noise  level  (5%  to  20%)  and  added  to  the 
noiseless  simulated  phase  of  Figure  2.  As  an  example,  for  a 
1%  noise  level,  an  input  noise  image  v  is  scaled  such  that 
=  0.01,  where  std  indicates  the  standard  deviation. 
Figure  2  (B)  shows  an  example  2D  cross  section  of  the  noisy 
simulated  phase  for  a  20%  noise  level. 

For  the  two  inverse  algorithms  an  optimal  value  of  the 
Tikhonov  factor  was  found  empirically,  and  used  for  each  of 
the  different  noise  levels.  A  value  of  lelO-5  was  used  for 
the  Quadratic  Regularization  reconstruction,  and  a  value  of 
5el0-6  for  the  TV  reconstruction.  The  effect  of  the  choice  of 
the  Tikhonov  value  is  discussed  later  in  this  Section. 

The  direct  algorithm  successfully  reconstructs  the  phase 
information,  showing  a  higher  conductivity  on  the  left  side 
of  the  domain,  and  a  vertical  transition  between  the  more 
conductive  and  less  conductive  regions.  The  gray,  3 -pixel  wide, 
band  present  at  the  boundary  of  the  image  is  an  artifact  of 
using  derivative  filters,  which  cannot  operate  in  proximity  of 
the  boundary;  the  derivative  filters  used  here  require  three 
pixels  per  side  of  the  considered  pixel.  The  inverse  quadratic 
algorithm  successfully  reconstructs  the  conductivity  profile 
showing  a  more  conductive  and  a  less  conductive  region  on 
the  left  and  right  sides,  respectively.  This  algorithm  provides 
a  similarly  smooth  transition  at  the  conductivity  boundary  as 
compared  to  the  direct  algorithm,  but  visually  is  more  stable 
in  the  presence  of  noise.  Total  Variation-based  reconstructions 
exhibit  a  sharp  transition  in  the  conductivity  distribution, 
resulting  in  a  better  overall  estimation  of  the  original  data.  This 
algorithm  is  also  the  least  sensitive  to  noise  for  this  dataset. 
TV  is  an  appropriate  image  prior  for  distributions  with  sharp 
transitions  like  the  one  used  for  these  tests,  and  therefore  likely 
to  produce  better  results  compared  to  other  methods.  All  the 
images  produced  by  the  inverse  approach  present  a  border  of 
one  pixel.  Pixels  on  the  boundary  are  not  estimated  as  they 
can  be  affected  by  the  boundary  condition  <p(r)  =  </>(r)meas. 

Figure  4  shows  the  effect  of  the  Tikhonov  factor  on  the 
reconstructed  images.  This  figure  shows  2D  cross-sections  of 
three-dimensional  reconstructions  with  the  Quadratic  Regu¬ 
larization  inverse  algorithm  for  a  noise  level  of  10%  and 
for  different  values  of  the  Tikhonov  factor,  which  are,  from 
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Fig.  1.  Simulated  domain  used  for  the  numerical  experiment:  a  20x20x20  mm  cube  with  a  conductivity  of  1  Sm-1  on  one  side  and  of  2  Sm-1  on  the 
other  was  generated  in  MATLAB.  The  simulation  of  the  MR  phase  within  the  cube  was  computed  for  a  3T  magnet  using  a  Finite  Difference  discretization 
with  1mm  resolution,  using  the  modeling  approach  developed  in  Section  II. 


Fig.  2.  2D  cross  section  in  the  horizontal  plane  of  the  3D  computed  MRI  phase  for  the  conductivity  distribution  shown  in  Figure  1.  Subfigure  (A)  shows 
the  noiseless  phase,  and  subfigure  (B)  the  noisiest  phase  (20%  noise  level)  fed  into  the  numerical  simulations.  As  expected  the  computed  phase  has  a  higher 
curvature  corresponding  to  the  more  conductive  region  (left  of  the  picture),  as  the  Laplacian  takes  larger  values  for  high  conductivity  regions  compared  to 
less  conducting  regions  (right  of  the  picture).  The  units  used  for  plotting  the  phase  are  radians,  and  the  range  is  -0.08  to  0.00  radians. 


left  to  right,  lelO-7,  lelO-6,  lelO-4,  and  lelO-3.  These 
values  bracket  the  optimal  value  of  lelO-5  used  in  the 
reconstructions  of  Figure  3.  The  two  left  reconstructions  are 
under-regularized  (i.e.  a  too  small  Tikhonov  value  results  in 
noisy  images).  The  two  right  reconstructions  are  instead  over¬ 
regularized  (i.e  a  too  large  value  of  Tikhonov  factor  results  in 
reconstructions  that  are  not  as  sensitive  to  noise,  but  in  which 
the  spatial  transitions  have  been  smoothed  excessively).  The 
value  of  lelO-5  used  in  Figure  4  therefore  represents  a  good 
compromise  between  sensitivity  to  noise  and  resolution. 

Conductivity  profiles  across  the  middle  of  the  domain  and 
passing  through  the  transition  point  were  extracted  from  the 
three  reconstructions  with  the  three  different  algorithms  for  the 
10%  noise  level  case  ,  and  are  shown  in  Figure  5).  The  orig¬ 
inal  conductivity  profile  (Figure  5 A)  used  for  generating  the 
synthetic  phase  data  shows  the  true  transition  each  algorithm 
attempts  to  recover.  The  bold  lines  depict  the  reconstructed 
profile,  while  the  true  conductivity  profile  is  shown  as  a  dotted 
line  in  Figures  5 A,  5B,  and  5C.  While  the  direct  approach 
(Figure  5 A)  and  the  inverse  approach  with  Quadratic  Regular¬ 
ization  (Figure  5B)  are  substantially  equivalent  in  terms  of  how 
fast  they  describe  the  transition  from  high  conductivity  to  low 
conductivity,  the  reconstruction  with  TV  regularization  (Figure 
5C)  shows  a  much  steeper  transition  in  the  reconstructed 
profile,  and  is  therefore  a  more  accurate  reconstruction  of  the 
step  conductivity  profile.  The  inverse  formulation  developed 
here  enables  one  to  select  regularization  functionals  that  are 


appropriate  for  the  problem  at  hand.  Besides  TV  and  similar 
edge-preserving  techniques  it  is  possible  to  envision  using  ad- 
hoc  functionals  that  incorporate  prior  structural  information 
extracted  from  other  anatomic  MR  images  [18],  [19]. 

VI.  Computational  Burden  Considerations 

While  the  developed  approach  offers  flexibility  in  terms  of 
the  regularization  functionals  available  and  does  not  require 
differentiating  the  noisy  input  images,  it  does  require  a  higher 
computational  load  compared  to  methods  based  on  direct 
differentiation.  This  computational  burden  can  be  measured 
in  terms  of  number  of  computations  and  amount  of  working 
memory  required.  We  will  show  that  the  computational  burden 
can  be  significant,  for  images  with  a  moderate  number  of 
pixels,  but  that  it  is  possible  to  reduce  it  by  subdividing 
the  image  domain  and  by  reconstructing  smaller  subdomains 
individually.  We  will  use  this  approach  for  reconstructing  the 
experimental  MRI  data  presented  in  Section  VII. 

Image  reconstruction  with  the  developed  inverse  formula¬ 
tion  requires  solving  the  forward  problem  (5),  using  for  exam¬ 
ple  a  finite-difference  scheme,  and  updating  the  conductivity 
or  permittivity  values  with  a  Gauss  Newton  update  as  in  (11). 
For  moderate  to  relatively  large  imaging  volumes  solving  (5) 
or  (11)  presents  challenges.  For  example,  an  imaging  volume 
of  256  x  256  x  256  pixels  requires  the  forward  problem 
to  be  solved  for  2563,  or  more  than  16  million  unknowns. 
Problems  of  this  size  can  be  efficiently  solved  on  a  desktop 
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Fig.  3.  2D  cross-sections  of  three-dimensional  reconstructions  of  synthetic  phase  data  with  noise  levels  of  5%,  10%,  15%,  and  20%,  associated  with  columns  1 
to  4,  respectively.  Each  row  in  the  figure  represents  a  reconstruction  with  a  different  algorithm:  the  top  row  shows  (for  different  levels  of  noise)  reconstructions 
with  our  own  implementation  of  the  algorithm  proposed  in  [3].  We  use  an  integration  volume  of  5  x  5  x  5  pixels,  and  estimate  derivatives  on  3  pixels 
per  side  of  the  central  pixel,  this  results  in  a  band,  3-pixels  wide,  around  the  image  that  cannot  be  reconstructed.  The  reconstruction  successfully  identifies 
the  conductivity  distribution,  showing  a  left-to-right  transition.  The  second  and  third  row  of  the  figure  represent  reconstructions  with  the  developed  inverse 
approach,  using  Quadratic  Regularization  (QR)  and  Total  Variation  (TV)  regularization,  respectively.  The  QR  algorithm  is  able  to  identify  the  left-to-right 
change  in  conductivity  and  to  describe  it  with  a  certain  degree  of  smoothness,  which  is  characteristic  of  this  type  of  regularization.  The  TV  algorithm  is  able 
to  identify  and  describe  with  pronounced  sharpness  the  left-to-right  change  in  conductivity.  In  this  specific  case  the  TV  algorithm  also  appears  particularly 
robust  to  noise.  For  both  QR  and  TV  algorithms  an  optimal  Tikhonov  factor  was  found  empirically,  and  maintained  constant  across  the  different  levels  of 
noise.  Specifically  a  value  of  lelO-5  and  of  5el0-6  was  used  respectively  for  the  two  algorithms.  With  the  parameters  used,  the  inverse  based  algorithms 


seem  to  fare  batter  in  the  presence  of  noise  compared  to  the  direct  algorithm, 
input  phase  data.  All  figure  are  represented  on  a  grayscale  ranging  from  0  to  2 


TABLE  I 

Computational  Time  For  Different  Subdomain  Configurations 


Config. 

Time  QR  Recon  [s] 

Time  TV  Recon  [s] 

Jac.  Mem.  [MB] 

2x1 

50.9 

684.8 

162 

2x2 

8.8 

232.9 

39 

3x2 

3.5 

133.3 

16 

3x3 

1.5 

85.9 
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computer  in  a  few  minutes  of  computation  time  using  special¬ 
ized  algorithms,  such  as  those  based  on  Algebraic  Multigrid 
Methods  [28].  These  algorithms  are  not  readily  available  in 
computing  environments  like  MATLAB,  and  require  the  user 
to  procure  specialized  libraries.  In  addition  to  the  time  required 
for  forward  solving,  there  is  a  particular  challenge  posed  by 
the  memory  utilization  required  to  solve  (11).  The  number 
of  rows  and  columns  in  the  Jacobian  matrix  J  is  equal  to 
the  number  of  unknowns  in  the  problem.  In  the  particular 
case  of  a  256  x  256  x  256  imaging  domain,  the  Jacobian 
would  be  a  2563  x  2563  matrix,  requiring  memory  allocation 
beyond  the  limits  of  any  personal  computer.  The  approach  we 
have  implemented  subdivides  the  imaging  domain  into  several 
smaller  subdomains  over  which  the  Jacobian  is  formed;  this 


as  it  would  be  expected  from  the  fact  that  they  do  not  need  to  differentiate 
Sm"1. 

significantly  reduces  the  computational  burden  and  memory 
requirements.  For  a  TV  x  TV  x  TV  imaging  domain,  the  size 
of  the  forward  problem  is  TV 3  and  the  memory  required  for 
storing  the  Jacobian  is  proportional  to  TV3* TV 3  =  TV6;  dividing 
the  domain  into  a  small  number  of  subdomains  immediately 
reduces  the  computational  burden. 

If  the  domain  is  subdivided,  for  example  into  cubic  blocks 
smaller  than  the  full  domain,  the  boundary  condition  (7)  can  be 
applied  at  the  interface  between  the  different  blocks.  The  for¬ 
ward  problem  can  then  be  solved  (5)  within  each  subdomain. 
Since  (7)  can  be  applied  everywhere  using  the  measured  phase 
values,  the  forward  problem  can  be  solved  on  a  subdomain 
independently  from  the  neighboring  subdomains,  and  image 
reconstruction  carried  out  on  each  subdomain  independently 
from  others. 

The  only  potential  dependence  between  neighboring  pixels 
is  introduced  at  the  interface  between  two  imaging  subdomains 
if  the  regularization  matrix  L  in  (11),  or  in  the  equiva¬ 
lent  Total  Variation  formulation  (17),  correlates  neighbouring 
pixels  across  the  two  different  subdomains.  In  our  current 
implementation,  we  do  not  regularize  pixels  across  different 
imaging  subdomains;  this  enables  us  to  run  fully  independent 
reconstructions  on  portions  of  the  full  domain.  For  real  MRI 
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Fig.  4.  Effect  of  the  Tikhonov  factor  on  reconstruction.  The  four  figures  show  2D  cross-sections  of  three-dimensional  reconstructions  of  synthetic  phase  data 
with  a  fixed  noise  level  of  10%  and  a  varying  Tikhonov  factor  of  a  value  of  lelO-7,  lelO-6,  lelO-4,  and  lelO-3,  from  left  to  right.  The  chosen  values 
bracket  the  optimal  value  of  lelO-5,  which  was  found  empirically  and  used  in  the  reconstructions  in  Figure  3.  The  values  lelO-7  and  lelO-6  result  in  an 
under-regularized  image  which  is  sensitive  to  noise.  The  values  of  lelO-4,  lelO-3  result  in  images  that  are  over-regularized,  where  sensitivity  to  noise  has 
been  greatly  reduced,  but  spatial  transitions  have  been  overly  smoothed.  All  figures  are  represented  on  a  grayscale  ranging  from  0  to  2  Sm- 1 . 


True  Conductivity 


1st  Derivative  Reconstruction 


Quadratic  Inverse  Reconstruction 


Total  Variation  Inverse  Reconstruction 


Fig.  5.  Subfigure  A  shows  a  plot  of  the  conductivity  value  on  a  horizontal  line  crossing  the  cube  of  Figure  1,  where  the  vertical  axis  represents  the 
conductivity  value  and  the  horizontal  axis  the  position  along  the  left-to-right  direction  inside  the  conductivity  cube.  Subfigures  B,  C,  and  D  represent  a  similar 
plot  for  the  reconstructed  conductivity  values,  for  a  noise  level  of  10%,  respectively  for  our  own  implementation  of  the  algorithm  proposed  in  [3],  for  the 
inverse  reconstruction  with  quadratic  regularization,  and  for  the  inverse  reconstruction  with  Total  Variation  regularization.  In  these  three  subfigures  the  bold 
continuous  line  represents  the  reconstructed  values,  and  the  thin  dotted  line  the  true  conductivity  value,  as  in  subfigure  A.  The  algorithm  in  subfigure  B  and 
C  have  a  similar  performance  in  terms  of  how  steeply  they  can  describe  the  sharp  transition,  while  the  algorithm  in  subfigure  D  (Total  Variation)  is  able  to 
describe  the  step  conductivity  change  much  more  accurately. 


data  (see  Section  VII),  we  split  the  imaging  domain  using 
this  sub-domain  technique  to  reduce  the  computational  burden 
associated  with  the  large  number  of  pixels  within  the  MR 
image  stack.  The  number  of  pixels  recorded  during  a  standard 
MR-EPT  scan  are  significantly  more  than  those  used  in  our 
numerical  simulation  (Section  V). 

To  substantiate  the  above  discussion,  we  report  the  compu¬ 
tational  time  and  the  required  amount  of  memory  for  storing 
the  Jacobian  matrix  for  different  subdomain  sizes  used  for 
reconstructing  the  phantom  in  Figure  9.  The  original  MRI 
image  (acquisition  details  described  in  Section  VII)  has  a  size 
of  144  x  144  x  144  pixels.  Within  this  volume  a  region  of 
interest  (ROI)  of  61  x  55  x  5  pixels  was  chosen  to  capture  the 
curved  detail  of  the  phantom  and  to  exclude  the  boundaries  of 
the  plastic  container  used  to  house  it.  A  magnitude  image  of 


the  selected  ROI  is  shown  in  Figure  9(B).  The  slices  of  the  ROI 
were  split  in  1,  2,  or  3  parts  along  the  vertical  and  horizontal 
directions,  giving  raise  to  a  number  of  smaller  subdomains 
for  image  reconstruction.  Precisely  we  used  the  following 
configurations:  2x1,  2x3,  3x2,  and  3x3,  where  the  first 
number  indicates  how  many  subdivisions  were  used  along  the 
vertical  direction  and  the  second  number  indicates  the  number 
of  subdivisions  along  the  horizontal  direction.  For  example, 
the  configuration  3x3  results  in  9  subdomains  in  total,  where 
slices  are  subdivided  with  a  3  by  3  grid  (see  Figure  10).  All 
the  information  (5  slices)  was  used  in  the  third  dimension. 

Table  I  reports  timing  and  memory  usage  information. 
All  computation  where  performed  on  a  workstation  with  a 
Intel  Xeon  3503  CPU  running  at  2.40GHz,  with  8GB  of 
memory,  using  the  Windows  7  Ultimate  -  64bit  operating 


9 


system;  algorithms  were  implemented  and  run  in  the  MAT- 
LAB  environment.  An  important  difference  exist  between  the 
algorithm  implementing  Quadratic  Regularization  and  and  the 
algorithm  implementing  Total  Variation  regularization  (TV): 
the  quadratic  update  equation  (11)  is  linear  with  respect  to 
<7  (the  Jacobian  does  not  depend  on  a ),  and  therefore  the 
term  [JT J  +  olLt L\~x  only  needs  to  be  computed  once, 
and  can  be  applied  later  to  any  subdomain  of  the  image  - 
a  fast  matrix-vector  multiplication.  The  Total  Variation  regu¬ 
larization  algorithm  instead  requires  an  iterative  cycle  on  each 
subdomain  of  the  image  (we  use  a  fixed  number  of  iterations, 
10  cycles  per  subdomain),  resulting  in  longer  reconstruction 
times  compared  to  Quadratic  Regularization  (QR).  For  both 
algorithms  computing  was  performed  in  double  precision, 
timing  and  memory  usage  information  is  therefor  relative  to 
an  8-byte  representation  of  floating  point  values.  The  coarser 
subdivision,  2x1,  requires  61  seconds  to  compute  with  QR 
and  685  seconds  with  TV,  using  162  Megabytes.  These  timing 
and  memory  requirements  reduce  respectively  to  1.5  and  86 
seconds  for  the  smaller  3x3  subdomain  configuration,  and 
to  7MB  of  memory  for  Jacobian  storage.  While  the  QR 
algorithm,  for  smaller  subdomains,  presents  a  computational 
burden  that  is  not  dissimilar  to  that  of  direct  approaches, 
TV  reconstruction,  as  implemented  in  our  algorithm,  still 
presents  a  burden  that  is  significant.  Use  of  TV  regularized 
reconstruction  is  therefore  a  tradeoff  between  the  benefits 
offered  by  this  type  of  functional  and  the  time  consumed  in 
the  reconstruction. 

VII.  Physical  Experiments 

The  proposed  MR-EPT  reconstruction  approach  is  demon¬ 
strated  on  two  phantom  datasets  acquired  using  a  Philips 
Achieva  3T  MR  scanner  using  a  Philips  SENSE  Flex  S  trans¬ 
mit/receive  coil.  A  first  conductivity  phantom  was  prepared 
by  slicing  a  12cm  x  7cm  x  5cm  block  of  gelatin  (a  = 
1.8  Sm-1)  into  slabs  of  different  thicknesses  (see  Figure  6). 
Specifically,  slice  thickness  of  20mm,  15mm,  10mm,  and  5mm 
were  formed  and  placed  into  a  polymer  housing.  The  slices 
were  positioned  inside  the  housing  and  adhered  to  the  walls  by 
slightly  heating  the  container;  the  heating  causes  a  thin  layer  of 
gelatin  to  melt  and  re- solidify  providing  adhesion.  Gelatin  was 
prepared  by  using  10%  by  weight  of  dry  porcine  gelatin,  and 
approximately  1%  by  weight  of  NaCl  to  adjust  conductivity 
to  the  desired  level.  A  saline  solution  with  a  conductivity  of 
4.1  Sm-1  was  produced  using  approximately  2%  NaCl  in 
weight  and  used  to  fill  the  gaps  between  the  slices  of  gelatin 
(not  shown  in  Figure  6),  to  provide  an  inclusion  (gelatin)  to 
background  (saline)  contrast.  In  addition,  approximately  1% 
by  weigh  of  copper  sulfate  (CuS04)  was  added  to  the  gelatin, 
but  not  to  the  saline,  in  order  to  provide  MR  contrast. 

Standard  MR  images  (e.g.  T1 -weighted  and  T2-weighted) 
were  not  expected  to  detect  the  electrical  properties  contrast 
between  the  gelatin  and  the  saline  solutions;  the  CuS04 
provides  MR  contrast  so  that  the  structure  of  the  phantom  can 
be  appreciated  in  standard  MR  images.  MRI  magnitude  images 
depict  a  high  intensity  (white)  where  the  gelatin  slabs  are 
located  (corresponding  to  the  high  concentration  of  CuS04), 


while  the  saline  solution  appears  as  a  low  intensity  (black) 
region  (see  Figure  6).  This  magnitude  MR  image  provides 
a  high  resolution  image  of  the  experimental  configuration 
to  which  reconstructed  MR-EPT  images  can  be  compared. 
A  standard  Spin  Echo  (SE)  MRI  sequence  was  utilized  for 
acquiring  the  magnitude  and  phase  data,  with  the  following 
settings:  resolution  of  80  x  80  x  20  pixels,  field  of  view 
(FOV)  of  160  x  160  x  60mm,  repetition  time  (TR)  =  600 
ms,  and  echo  time  (TE)  =  7.79  ms.  Four  temporal  averages 
were  used  in  order  to  improve  the  SNR  for  a  total  acquisition 
time  of  389  seconds.  MR  phase  information  was  captured 
together  with  amplitude  information  and  used  for  reconstruct¬ 
ing  conductivity  images  with  the  inverse  implementation  and 
with  the  direct  method  proposed  by  Katscher  et.  al.  [3]. 
For  the  inverse  formulation,  both  quadratic  regularization 
and  TV-based  regularization  approaches  were  used.  For  all 
reconstructions,  data  extracted  from  a  central  portion  of  the 
phantom,  consisting  of  a  58  x  22  x  7  pixel  volume,  was 
used.  All  three  reconstruction  approaches  exhibit  regions  of 
high  and  low  conductivity  corresponding  to  the  saline  and 
gelatin,  respectively  (see  Figure  7).  The  same  grayscale  was 
selected  for  each  of  the  different  algorithms,  ranging  from 
0  Sm-1  (black)  to  4  Sm-1  (white).  It  is  important  to  note 
that  gelatin  slices  present  a  lower  conductivity  (1.8  Sm-1) 
compared  to  the  saline  solution  (4.1  Sm-1),  making  the  gelatin 
slices  appear  in  a  darker  gray  level.  This  is  in  contrast  to 
MR  magnitude  images  (Figure  6),  where  gelatin  slices  appear 
brighter  due  to  the  presence  of  (CuS04).  These  contrast 
differences  observed  between  MR  magnitude  images  (Figure 
6)  and  MR-EPT  reconstructions  (Figure  7)  help  to  validate 
MR-EPT ’s  dependence  on  conductivity  and  not  on  typical  MR 
contrast  mechanisms  (e.g.  CuS04). 

All  three  reconstruction  approaches  correctly  identify  the  lo¬ 
cation  of  the  gelatin  strips  within  the  phantom  (Figure  7).  The 
reconstruction  using  our  implementation  of  the  direct  method 
developed  by  Katscher  et.  al.  [3]  has  a  border  of  3  pixels 
over  which  the  conductivity  values  are  not  reconstructed.  This 
artifact  stems  from  using  Savitzky-Golay  filters  with  3  pixels 
on  each  side  of  the  central  pixel  of  interest  to  estimate  the 
derivatives  at  these  locations.  In  the  inverse  formulation,  a 
single  pixel  around  the  border  is  used  to  match  the  measured 
phase  with  the  Dirichlet  boundary  condition  (7)  and  is  not 
reconstructed;  all  other  interior  pixels  are  fit  to  the  data.  In  this 
particular  experimental  configuration,  the  reconstruction  with 
quadratic  regularization  (Figure  7B)  provides  slightly  more 
contrast  as  compared  to  the  reconstruction  based  on  integration 
of  first  derivatives  (Figure  7A).  Much  more  noticeable  is 
the  difference  between  the  reconstruction  based  on  Total 
Variation  (TV)  (Figure  1C)  regularization  and  the  other  two 
reconstruction  approaches.  TV  represents  a  good  prior  for 
describing  the  sharp  conductivity  transitions  that  are  present  at 
the  gelatin- saline  interface;  using  this  form  of  regularization 
significantly  enhances  the  reconstructed  images.  The  recon¬ 
structed  conductivity  (Figure  1C)  geometrically  aligns  with 
the  phantom  configuration  (as  viewed  visually  and  in  MR 
magnitude  images)  in  that  the  high  conductivity  gelatin  regions 
(dark  stripes)  correspond  to  high  intensity  MR  magnitude 
regions  (light  stripes).  The  TV-based  regularization  accurately 
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Fig.  6.  Top:  photographic  representation  of  the  conductivity  phantom.  A  conductivity  phantom  was  generated  by  slicing  a  slab  of  gelatin  in  slices  of  20mm, 
15mm,  10mm,  and  5mm  of  thickness.  The  slices  were  positioned  inside  a  polymer  housing  and  secured  to  it  by  slightly  heating  the  housing  and  letting  the 
gelatin  in  contact  with  the  housing  walls  slightly  melt  before  re-solidification.  The  phantom  was  then  filled  with  a  saline  solution  (not  shown).  The  gelatin 
had  a  conductivity  of  1.8  Sm-1  and  the  saline  solution  of  4.1  Sm-1.  Copper  sulfate  (CuS04)  was  added  to  the  gelatin,  so  that  the  gelatin  slabs  would  show 
with  a  different  intensity  in  traditional  MR  amplitude  images.  Bottom:  an  MR  amplitude  image  shows  the  phantom.  Gelatin  slices  appear  in  a  brighter  gray 
due  to  the  inclusion  of  Copper  Sulfate  (note  that  gelatin  slices  have  a  lower  conductivity  value  compared  to  the  saline  solution  (1.8  Sm-1  versus  4.1  Sm-1), 
and  therefore  they  appear  darker  in  the  MR-EPT  reconstructions  of  Figure  7). 
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Fig.  7.  MR-EPT  Reconstructions  of  MR  data  acquired  on  the  phantom  of  figure  6.  Subfigure  A  shows  a  “traditional”  MR-EPT  reconstruction  with  our 
implementation  of  the  algorithm  proposed  in  [3].  A  border  of  3  pixels  is  present  all  around  the  boundary  of  the  image,  as  derivatives  are  estimated  using 
three  pixel  values  per  each  side  of  the  current  pixel,  therefore  reconstruction  can  occur  only  for  pixels  that  are  at  a  distance  of  3  pixels  from  the  domain 
boundary.  The  algorithm  reconstructs  correctly  the  alternating  conductivity  values,  corresponding  to  the  high  conducting  saline  solution,  and  less  conducting 
gelatin  slices.  Gelatin  slices  appear  in  fact  darker,  as  opposed  to  the  lighter  gray  in  which  they  appear  in  the  MR  magnitude  image  in  Figure  6B,  as  the 
contrast  mechanism  in  MR-EPT  is  based  on  the  electrical  properties,  while  the  MR  amplitude  image  is  sensitive  to  the  presence  of  Copper  Sulfate.  Subfigure 
B  shows  a  MR-EPT  reconstruction  using  the  developed  inverse  formulation  with  quadratic  regularization,  and  Subfigure  C  a  MR-EPT  reconstruction  using 
the  developed  inverse  formulation  with  Total  Variation  regularization.  All  figure  are  represented  on  a  grayscale  ranging  from  0  to  4  Sm_  1 . 
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Phase  Fitting  Process  -  Quadratic  Algorithm 


Fig.  8.  Phase  fitting  process.  This  figure  shows  how  the  MR-EPT  inverse  reconstruction  algorithm  is  fitting  the  phase  data.  The  plot  on  the  left  shows  results 
relative  to  the  QR  algorithm,  and  the  plot  on  the  left  result  relative  to  the  TV  algorithm.  The  dash-dotted  line  at  the  top  of  the  plots  represents  the  computed 
phase  value  across  a  numerical  phantom  of  Figure  6  corresponding  to  a  uniform  initial  distribution  of  conductivity.  The  solid  bold  line  at  the  bottom  of 
the  plots  represents  the  true  measured  phase  across  the  phantom.  The  dips  and  peaks  of  the  measured  phase  represent  phase  changes  corresponding  to  the 
alternating  high-low  conductivity  values  encountered  across  the  gelatin  stripes  and  saline  volumes  that  constitute  the  phantom  (see  Figure  6A).  The  dash-dotted 
line  at  the  bottom  of  the  plots  represents  the  fitted  phase,  resulting  from  the  inversion  procedure.  Despite  reconstructed  images  looking  quite  different  for 
the  QR  and  TV  algorithms,  the  fitted  phases  show  minor  differences.  This  is  to  be  expected,  as  different  values  of  reconstructed  conductivity  can  stem  from 
minor  curvature  changes  in  the  phase.  For  both  algorithms  the  fitted  phase  follows  closely  the  general  trends  of  the  measured  phase;  regularization  techniques 
used  in  the  algorithm  help  to  ensure  that  the  model  does  not  “overfit”  the  small  scale  phase  perturbations  that  represent  noise  in  the  data  (e.g.  small  dip  near 
pixel  position  20). 


Fig.  9.  Curved  phantom  photographic  image  (left)  and  MRI  magnitude  image  (right).  This  phantom  was  build  with  a  similar  method  to  the  phantom  shown 
in  Figure  6,  by  creating  a  gelatin  slab  which  is  immersed  in  a  saline  solution.  The  gelatin  slab  has  a  conductivity  of  of  1.8  Sm-1  and  Copper  sulfate  (CuS04) 
was  added  to  it,  to  generate  contract  in  the  MRI  magnitude  image.  The  saline  solution  has  a  conductivity  of  4.1  Sm-1.  In  part  (A)  of  the  figure  is  indicated 
a  white  edge  that  approximates  the  perimeter  of  a  region  of  interest  which  was  selected  for  image  reconstruction.  Part  (B)  of  the  figure  shows  the  MRI 
magnitude  corresponding  to  the  region  of  interest  used  in  the  reconstructions  of  Figure  10. 


identifies  the  steep  transitions  that  are  characteristic  of  the 
phantom,  while  the  direct  and  quadratic-based  inverse  formu¬ 
lation  do  not  recover  the  steep  transitions.  Both  quadratic  and 
TV-regularized  inverse  approaches  to  reconstruction  (Figure 
7  B  and  C)  exhibit  a  conductivity  discontinuity  through  the 
center  of  the  image.  This  discontinuity  arises  from  having 
subdivided  the  reconstruction  domain  into  an  upper  and  lower 
subdomain  for  the  purpose  of  speeding-up  reconstruction  as 
discussed  in  Section  VI.  In  the  current  implementation  of  this 
approach  the  different  subdomains  in  which  an  image  can  be 
divided  are  reconstructed  in  a  completely  independent  manner 
(i.e.  we  do  not  introduce  correlation  between  neighbouring 
pixels  at  the  interface  of  the  subdomains  in  the  regularization 
matrices),  and  results  in  small  discontinuities  observed  in  the 
images.  Using  a  scheme  that  incorporates  subdomain  overlap 


as  presented  in  [29]  will  be  explored  in  the  future  to  reduce 
or  eliminate  this  artifact. 

The  fitting  process  resulting  from  the  inverse  approach 
to  reconstruction  can  be  evaluated  by  observing  how  the 
fitted  phase  progresses  from  the  initial  guess  to  the  final 
reconstructed  phase  (Figure  8).  In  this  case,  the  estimated 
phase  (y-axis)  is  plotted  as  a  function  of  distance  (x-axis)  along 
a  trajectory  passing  horizontally  through  the  striped  phantom 
for  both  the  QR  and  TV  algorithms.  The  dips  and  peaks  in 
the  phase  data  are  caused  by  the  varying  conductivity  values 
across  the  phantom  gelatin  slabs  and  saline  volumes.  True 
measured  phase  is  indicated  with  a  solid  line,  and  the  dash- 
dotted  line  at  the  top  of  the  plots  represents  the  computed 
phase  for  an  initial  uniform  conductivity  distribution  of  1 
Sm-1,  while  the  dash-dotted  line  at  the  bottom  of  the  plots 
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Fig.  10.  2D  cross-sections  of  three-dimensional  reconstructions  of  the  phantom  of  Figure  10.  This  figure  demonstrates  the  ability  of  the  inverse  algorithms 
to  reconstruct  curved  boundaries,  and  also  serves  to  demonstrate  the  possibility  of  speeding  up  image  reconstruction  by  splitting  the  image  domain  in  smaller 
subdomains  to  be  reconstructed  independently.  The  left  column  represents  reconstructions  with  the  inverse  algorithm  with  Quadratic  Regularization  and  the 
right  column  reconstructions  using  Total  Variation  regularization.  Each  row  represents  a  different  subdomain  splitting  for  image  reconstruction,  as  discussed 
in  Section  VI.  In  the  top  row  MRI  slices  have  been  split  in  two  parts  vertically,  indicated  as  2  x  1  splitting,  in  the  second  row  the  slices  have  been  split  in  two 
parts  vertically  and  horizontally,  indicated  as  2  x  2  splitting,  and  the  third  and  fourth  represent  respectively  3x2  and  3x3  splittings.  Splitting  the  original 
image  domain  in  smaller  parts  allows  for  a  faster  reconstruction  as  discussed  in  Section  VI  and  reported  in  Table  I.  Some  discontinuities  are  present  at  those 
image  locations  where  two  different  subdomains  meet.  This  results  from  the  fact  that  different  subdomains  are  treated  as  separate  by  the  reconstruction  and 
no  continuity  is  enforced.  In  future  work  we  intend  to  introduce  some  correlation  between  neighbouring  subdomains,  through  regularization  techniques,  which 
should  reduce  or  eliminate  these  discontinuities.  In  both  sets  of  reconstructed  images,  though  some  minor  discontinuities  are  present  at  the  interface  between 
subdomains  both  and  inverse  QR  and  inverse  TV  algorithms  are  able  to  correctly  reproduce  the  curved  interface.  All  figures  are  represented  on  a  grayscale 
ranging  from  0  to  4  Sm-1. 
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represents  the  computed  phase  after  the  fitting  is  complete. 
The  computed  phase  closely  tracks  the  measured  phase,  with 
the  exception  of  small-scale  features  associated  with  noisy 
measurements.  The  regularization  helps  to  enforce  an  accurate 
fit  of  the  measured  phase,  but  to  avoid  fitting  small  scale 
variations  that  are  typically  associated  with  measurement  noise 
and  errors. 

A  second  conductivity  phantom  was  prepared  using  the 
same  method  described  above  for  the  striped  phantom  of 
Figure  6.  In  this  case  though  a  portion  of  a  round  slab  of 
gelatin  was  cut  and  placed  in  a  corner  of  a  polymer  housing  as 
shown  in  Figure  9.  This  second  physical  phantom  was  used  to 
evaluate  the  inverse  reconstruction  algorithms  on  a  curvilinear 
geometry.  The  striped  phantom  of  Figure  6  is  well  suited  for 
reconstruction  with  TV  regularization  algorithms,  as  it  presents 
straight  step  changes,  aligned  with  one  axis  of  the  image,  that 
can  be  accurately  reconstructed  by  TV  algorithms.  The  curved 
boundary  of  the  phantom  in  Figure  9  poses  a  challenge  to  TV 
based  algorithms,  as  they  are  know  to  tend  to  reconstruct  round 
boundaries  with  a  staircase  appearance.  The  reconstructions  of 
Figure  10  demonstrate  that  both  the  QR  and  TV  algorithms 
are  able  to  capture  and  reproduce  the  curved  nature  of  the 
boundary  between  gelatin  and  water.  Both  algorithms  exhibit 
some  artifact  at  the  interface  between  subdomains,  since  we 
have  not  enforced  any  contiguity  between  the  reconstructed 
values  from  different  subdomains.  TV  regularization  seems  in 
this  case  to  suffer  less  from  the  staircase  effect  compared  to 
application  in  other  tomographic  applications.  For  example  in 
Electrical  Impedance  Tomography  the  staircase  effect  seems 
to  be  more  pronounced  [25].  We  believe  this  might  be  due  to 
the  fact  that  in  MR-EPT  data  is  measured  everywhere  in  the 
domain  and  not  only  at  the  boundary  as  in  other  techniques. 
This  might  help  to  drive  the  reconstruction  towards  more 
realistic  results. 

VIII.  Conclusions 

A  novel  approach  to  MR-EPT  image  reconstruction  based 
on  an  inverse  formulation  has  been  developed.  The  approach  is 
based  on  the  observation  that  the  central  equations  of  MR-EPT 
can  be  inverted  and  used  to  describe  a  forward  model,  which 
in  turn  can  be  used  in  an  inverse,  data  fitting  approach  to  MR- 
EPT  reconstruction.  This  approach  is  valid  for  reconstruction 
of  conductivity  from  measured  phase  information  and  for 
reconstruction  of  permittivity  from  measured  H+  amplitude 
information,  though  in  the  present  manuscript  we  develop 
and  demonstrate  this  approach  only  for  reconstruction  of 
conductivity.  A  forward  model  for  computing  phase  infor¬ 
mation  from  conductivity  data  is  presented  (Section  II).  In 
addition,  an  inverse  approach  based  on  quadratic  regularization 
and  Jacobian  computation  has  been  developed  for  solving 
the  MR-EPT  image  reconstruction  problem  (Section  IV- A). 
An  inverse  formulation  using  Total  Variation  as  a  functional 
for  regularization  has  been  developed  (Section  IV-B);  this 
approach  enables  the  reconstruction  of  sharper  conductivity 
transitions  within  the  imaging  domain.  Numerical  experiments 
were  conducted  to  validate  the  developed  inversion  approaches 
in  the  presence  of  synthetic  noisy  data  (Section  V).  A  method 


for  splitting  the  imaging  domain  into  subdomains  has  been 
developed  to  reduce  the  computational  burden  arising  from  the 
proposed  inverse  approach  (Section  VI);  this  implementation 
requires  less  memory  utilization  and  fewer  computations  lead¬ 
ing  to  faster  reconstructions.  Subdividing  the  imaging  domain 
into  smaller  subdomains  results  in  significant  performance  / 
memory  utilization  gains,  as  both  the  memory  and  computation 
time  are  linked  to  the  size  of  the  Jacobian  matrix  (which 
grows  with  the  6th  power  of  the  size  of  the  imaging  domain). 
Successful  reconstructions  were  obtained  on  laboratory  data 
collected  from  a  phantom  experiment  that  included  alternating 
regions  of  high  and  low  conductivity,  and  from  a  phantom 
with  a  curved  boundary  (Section  VII).  Both  the  inverse 
approach  using  quadratic  regularization  and  the  approach 
using  Total  Variation  regularization  accurately  identified  the 
position  of  the  gelatin  slabs;  TV  regularization  more  accurately 
reconstructed  the  steep  conductivity  transitions  present  at  the 
gelatin- saline  interface  of  the  phantoms  imaged.  Benefits  of 
this  reconstruction  method  include:  1)  no  differentiation  is 
required  on  the  input  data,  therefore  decreasing  sensitivity  to 
noise  compared  to  other  methods  described  in  the  literature; 
2)  different  regularization  functionals  can  be  implemented, 
depending  on  the  expected  distribution  of  the  parameters 
to  be  reconstructed.  This  capability  is  particularly  useful  in 
the  context  of  reconstructing  biomedical  data,  in  such  cases 
custom  regularization  functionals  can  be  constructed  ad-hoc 
to  incorporate  prior  anatomical  information  and  ultimately  en¬ 
hance  the  quality  and  robustness  of  the  reconstructed  images. 
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ABSTRACT 


Determining  prostate  cancer  grade  is  critical  for  determining  optimal  treatment  strategies  for  men 
diagnosed  with  prostate  cancer.  MR-based  approaches  have  been  explored  as  a  modality  to  aid  clinicians 
in  properly  staging  men  with  prostate  cancer.  Here  we  explore  the  feasibility  of  using  an  MR-based 
imaging  modality  capable  of  mapping  the  electrical  conductivity  of  tissue  for  prostate  cancer  imaging. 
This  imaging  modality,  termed  Magnetic  Resonance  -  Electrical  Properties  Tomography  (MR-EPT),  is 
based  on  computationally  manipulating  phase  data  recorded  during  an  MR  scan  to  reconstruct  an  image 
of  the  electrical  conductivity.  In  the  present  work,  the  MR-EPT  reconstruction  method  employed  is  based 
on  a  novel  inverse  problem  formulation  recently  developed.  The  method  offers  three  main  advantages 
over  a  direct  inversion  approaches  typically  used  in  MR-EPT:  a)  no  spatial  differentiation  is  needed, 
reducing  the  impact  noise  has  on  image  formation,  b)  a  regularization  term  enables  user  control  of  the 
resolution  of  the  reconstructed  data,  and  3)  a  prior  information  (e.g.  anatomic  MR  images)  can  be  used  to 
help  guide  the  reconstruction.  The  feasibility  of  using  MR-EPT  for  prostate  imaging  is  explored  through 
use  of  gelatin  and  saline  phantoms,  biological  tissue  phantoms,  and  scans  of  ex  vivo  prostate  tissue. 


1.  INTRODUCTION 


MR-EPT  (Magnetic  Resonance  Electrical  Properties  Tomography)  is  a  technique  for  imaging/estimating 
the  electrical  properties  of  tissues  based  on  MRI  technology  without  the  need  for  electrodes  as  are  needed 
for  Electrical  Impedance  Tomography  (EIT)  or  Magnetic  Resonance  EIT  (MR-EIT).  MR-EPT  was 
originally  proposed  by  Katscher  et  al.  [1,  2]  in  2009  and  a  recent  review  of  the  technique  has  been  offered 
by  Zhang  et  al.  in  [3]. 

MREPT  is  based  on  the  fact  that  the  RF  B1  field,  used  for  flipping  the  magnetic  moments  of  the  protons, 
interacts  with  the  body  to  be  imaged  in  the  scanner;  as  a  result  its  magnitude  and  phase  gets  altered.  With 
appropriate  MRI  sequences  this  altered  B 1  field  data  can  be  gathered  (in  space,  and  within  the  body  to  be 
imaged).  Based  on  this  data  it  is  possible  to  reconstruct  the  distribution  of  electrical  conductivity  and 
permittivity  of  the  body  being  imaged.  Obtaining  maps  of  conductivity  and  permittivity  has  diagnostic 
and  clinical  value. 

MR-EPT  applied  to  prostate  imaging  (cancer  data,  outcomes)  -  [references  Ryan  and  Andrea]. 

2.  METHODS 

2.1  Inverse  problem  formulation 

The  general  approach  for  conductivity  imaging  with  MR-EPT  is  to  obtain  a  phase  image  and/or  B1  map 
image  of  the  RF  field  produced  using  specific  pulse  sequences.  The  data  contained  in  the  image  is  utilized 
to  estimate  the  conductivity  distribution  of  the  regions  within  the  image.  A  number  of  methods  have  been 
proposed  for  this  manipulation:  one  approach  requires  second  derivatives  to  be  computed  from  the  phase 
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data,  but  the  method  is  prone  to  noise  amplification;  other  approaches  have  included  algorithms  that  lower 
this  requirement  to  first  derivatives,  thus  reducing  sensitivity  to  noise.  ADD  REFS! 

This  work  presents  an  alternative  method  that  solves  the  MR-EPT  problem  using  an  inverse  problem 
formulation  that  does  not  require  differentiating  the  input  image,  as  described  by  [4]. 

Paragraph  from  EIT  2014: 


*The  electrical  conductivity  a  can  be  shown  to  be  proportional  to  the  Laplacian  of  the  phase  of  the 
transmit  B1  field: 


a(r)  ~  ^A<KH+(r)).  (1) 

The  inverse  is  true  as  well:  if  o(r)  is  known,  the  phase  can  be  obtained  by  solving  Ac|)  =  oo[ia(r). 

Using  an  iterative  inverse  formulation  approach,  the  updated  value  of  a  is  given  by  anew  =  a  +  6a  where 

5a  =  (JTJ  +  aLTL)_1  (jT(c|)(H  +  (r))  -  <|>  +  aLTLa)).  (2) 

Here  J  is  the  Jacobian  of  the  conductivity  to  phase  mapping,  L  is  a  regularization  matrix,  and  a  is  a 
regularization  parameter  used  to  stabilize  the  inversion. 

We  have  implemented  this  inversion  using  two  different  regularization  terms: 

a)  a  quadratic/Laplacian  approach  and 

b)  a  Total  Variation  functional  approach  [5,6].  A  Primal  Dual  Interior  Point  Method  optimization  scheme 

is  used  for  the  Total  Variation  approach,  which  produces  images  with  sharper  contrasts  at  boundaries. 

* 

TO  ADD  HERE! 

2.2  Data  acquisition 


2.2.1.  Phantoms 

To  validate  the  proposed  approach  experimentally  a  number  of  phantoms  comprised  of:  gelatin  and 
saline,  gelatin-gelatin,  gelatin  and  a  playdough  inclusion,  have  been  developed  and  imaged.  The  gelatin 
used  for  the  phantoms  to  be  presented  below  is  porcine  gelatin  (Sigma-Aldrich,  300  bloom).  Saline 
solution  (deionized  sterile  water  and  table  salt  (NaCl)  in  various  proportions  has  been  used  as  a  phantom 
and  also  to  vary  the  conductivity  of  gelatin,  by  adding  it  in  the  process  of  making  the  gelatin.  When 
desired,  cupric  sulphate  (Sigma-Aldrich)  was  added  for  MR  contrast.  An  electrical  conductivity  meter 
(company?)  has  been  used  to  measure  the  electrical  conductivity  of  any  solution  used  in  the  process  of 
making  the  phantom  (i.e.  before  adding  gelatin,  for  the  gelatin  phantom).  A  detailed  description  of  the 
composition  and  geometry  of  the  phantoms  used  in  the  current  study  follows: 

a)  Gelatin  and  saline  phantom  (aka  resolution  phantom) 
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Figure  1.  Gelatin  and  saline  phantom:  geometry  description  and  phantom. 


A  block  gelatin  phantom  (10%  gelatin,  1%  NaCl,  -5 5x  125x75mm)  was  constructed  with  three  rows  of 
circular  wells  with  increasing  diameters  (5,  10,  15mm).  The  gelatin  block  was  cured  for  24  hours  at  4°C. 
Prior  to  imaging,  each  series  of  wells  was  filled  with  saline  solutions  with  increasing  conductivities  (3,5, 
8  S/m,  as  measured  with  a  conductivity  meter  in  solution)  (Figure  1).  To  visually  differentiate  easily 
between  the  saline  solutions  filling  each  row  of  wells,  within  the  magnitude  MR  images,  different  minute 
amounts  of  cupric  sulfate  were  added  to  the  solutions. 

b)  Gelatin  phantom  with  two  regions 

A  gelatin  phantom,  with  the  geometry  of  a  real  prostate  (-50mm  in  diameter)  -  a  silicone  mold  created 
from  an  actual  prostate  imaged  geometry  was  used.  The  phantom  was  comprised  of  two  regions  of 
different  conductivities,  as  seen  in  Figure  2. 


Figure  2.  Gelatin  phantom  with  two  regions 

The  mold  was  initially  filled  halfway  with  gelatin  of  one  conductivity  and  after  the  gelatin  set  and  cooled 
(~lh  at  4°C),  the  second  gelatin  solution  was  poured  to  fill  the  mold.  The  phantom  was  cured  for  24  hours 
at  4°C.  To  differentiate  visually  between  the  two  conductivity  regions,  food  coloring  was  added  to  each 
saline  solution.  The  solutions,  before  the  addition  of  the  gelatin,  had  a  measured  conductivity  of  1.6  S/m 
and  3.4  S/m,  respectively  for  the  pink/right  region  and  blue/left  region  in  Figure  2.  Prior  to  imaging,  the 
phantom  was  submerged  in  a  bath  of  deionized  water,  supported  and  kept  in  place  by  paper  rings  that 
provided  no  interference  to  the  MR  signal  and  minimal  ringshaped  contact  to  the  phantom. 

c)  One  inclusion  gelatin  phantom 
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Figure  3.  One  inclusion  gelatin  phantom 


A  third  phantom  (Figure  3),  prostate-shaped,  with  one  5mm  inclusion  (made  of  play  dough,  for 
conductivity  contrast)  was  constructed  in  a  similar  manner  with  the  phantom  described  above.  The  saline, 
colored  with  food  coloring,  used  to  make  the  gelatin,  had  a  measured  conductivity  of  2 S/m.  The  phantom 
was  cured  for  24  hours  at  4°C.  Prior  to  imaging,  the  phantom  was  submerged  in  a  bath  of  deionized  water, 
supported  and  kept  in  place  by  paper  rings  in  a  similar  manner  with  what  is  shown  in  Figure  2. 

NOTE:  this  phantom  was  scanned  with  3mmx3mm(x3mm  for  the  voxel)  resolution.  The  previous  have  a 
2x2x3pixel  size!  It  was  in  the  interest  of  reducing  scan  time  in  preparation  for  the  prostate  scans. 

d)  Two  inclusions  gelatin  phantom 

Two  5  mm  diameter  inclusions  constructed  from  play  dough  to  provide  significant  conductivity  contrast, 
were  embedded  in  a  gelatin  prostate  with  a  conductivity  of  approximately  2  S/m,  in  a  similar  fashion  as 
described  for  phantom  c),  and  shown  in  Figure  3. 

NOTE:  voxel  is  2  mm  x  2  mm  x  3mm 

NEED  PHOTO? 


e)  Meat  phantom 

The  phantom  was  an  actual  block  of  pork  meat  (approximately  70x70x50mm),  as  seen  in  with  Figure  4, 
immersed  in  deionized  water  prior  to  being  imaged.  The  goal  of  reconstructing  this  phantom  was  to 
accurately  differentiate  between  the  muscle  fibers  and  fat  (a  layer  of  variable  thickness  on  the  top  region). 

NOTE:  this  phantom  was  scanned  with  3mmx3mm(x3mm  for  the  voxel)  resolution. 
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2.2.2.  Ex-vivo  prostate 


An  approved  IRB  protocol  for  imaging  excised  prostate  tissue  was  obtained  and  a  number  of  male 
subjects  consented  to  participate  in  this  study.  An  ex-vivo  prostate  was  imaged  within  an  hour  from 
excision,  in  order  to  be  ready  for  the  subsequent  pathology  analysis.  Upon  receiving  it,  the  prostate  was 
placed  in  an  acrylic  tub  filled  with  deionized  water  (low  conductivity)  and  supported  in  place  by  paper 
cylinders  (as  shown  in  Figure  2). 

Candidates:  PI 08  (big  tumor,  one  MREPT  image  that  looks  promising),  Pill,  P110. 

To  ask  Aditya:  age  of  the  specimen. 

Do  we  mention  the  type  of  cancer? 


2.2.3.  Imaging  Protocol 


QUESTIONS:  Use  the  SE  scans  for  all?  For  phantom  c)  -  we  have  a  TSE  scan,  but  have  data  for  similar 
phantoms  SE. 

Different  pixel  size? 

Do  I  describe  this  as  a  process  to  get  to  the  best  imaging  sequence,  as  it  was  historically  done,  or  do  I  use 
the  best  sequence  we  have?  Might  need  to  redo  some  phantoms. . . 

The  phase  data  necessary  to  reconstruct  the  conductivity  is  collected  on  Philips  Achieva  3T  scanner, 
using  a  Spin  Echo  (SE)  multi-slice  sequence  (TR=800,  TE=20)  with  a  flip  angle  of  90°. 

We  try  to  acquire  16  signal  averages  (NSA=16)  in  both  cases  in  the  interest  of  achieving  a  high  SNR. 
However  if  the  scan  duration  is  too  long  (>lh),  especially  in  the  TSE  sequence  case,  we  reduce  the  signal 
averages  (NSA=8).  The  data  is  acquired  with  a  flex  coil  (one  transmit,  one  receive  channel);  no 
endorectal  coil  is  used  clinically  at  DHMC.  In  all  scans,  for  the  ease  of  comparison  with  subsequent 
clinical  data,  we  use  a  patient  position  used  for  MR  scans  of  the  prostate:  head-first  supine,  with 
‘transverse’  slices  (from  seminal  vesicles  to  the  apex  of  the  prostate)  (fold-over  direction  AP,  fat  shift 
direction  L).  The  current  preferred  voxel  size  is  2x2x3mm  (a  3mm  imaging  slice  thickness  is  clinically 
used).  The  corresponding  reconstructed  voxel  size  is  1.63xl.63x3mm.  The  field  of  view  (FOV)  and 
number  of  slices  influence  the  scan  duration  as  well.  For  example  for  the  Gelatin  phantom  with  two 
regions  (phantom  b),  with  the  above  mentioned  parameters,  a  stack  of  25  slices  and  a  FOV  of 
130x130x75mm  lead  to  a  total  scan  duration  of  1:00:16  for  the  TSE  sequence  and  27:48.8  for  the  SE 
sequence.  Although  the  TSE  tends  to  result  in  sharper  contrasts  in  the  reconstructed  conductivity  maps, 
its  duration  makes  it  difficult  to  use  for  the  future  in  vivo  study,  so  the  SE  sequence  will  be  used  instead. 

3.  RESULTS 

Do  we  show  all  recons  or  only  3D  (laplacian  and  inverse,  NO  TV  3D?)? 
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4.  CONCLUSIONS 
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